clear all %#ok
close all
clc



disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGE 1 ESTIMATION:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
%%%%IMPORTANT NOTE: AAA_FILE01_ThetaEDriver.m includes a small 
%%%%correction to the data that is needed before we compute descriptive stats...
run AAA_FILE01_ThetaEDriver.m 

%%



AchievementData             = dataset ; %#ok
AchievementData.id          = g2.student_id ;
AchievementData.Q           = nan(length(AchievementData.id),1) ;  %%%%Total output
AchievementData.T           = nan(length(AchievementData.id),1) ;  %%%%This will contain all cumulative work time
AchievementData.ThetaE      = nan(length(AchievementData.id),1) ;  %%%%This will contain the subjects estimated ThetaE type (when available)
AchievementData.contract    = nan(length(AchievementData.id),1) ;  %%%%This will contain the numbering of the subject's wage contract (being 1, 2, or 3)
AchievementData.b           = nan(length(AchievementData.id),1) ;  %%%%This will contain the base wage for the subject's assigned wage contract
AchievementData.w           = nan(length(AchievementData.id),1) ;  %%%%This will contain the marginal wage for the subject's assigned wage contract
AchievementData.Tfirstunit  = zeros(length(AchievementData.id),1) ;  %%%%This will store work time on the first unit for active students, and all work time for marginal students
AchievementData.Tunfinished = zeros(length(AchievementData.id),1) ;  %%%%This will store work time on the first unit for active students, and all work time for marginal students

studlistFULL     = EPointEst.studlistFULL ; %%%%This is the student id list for everyone who completed at least 1 learning task 
studlistACTIVE   = EPointEst.studlist ;     %%%%This is the student id list for everyone who completed at least 2 learning tasks 
studlistMARGINAL = EPointEst.marginallist ; %%%%This is the student id list for everyone who completed exactly 1 learning task 
for ii=1:length(AchievementData.id) 
    if sum(studlistACTIVE==AchievementData.id(ii))+sum(studlistMARGINAL==AchievementData.id(ii))==0 
        %%%%This means the student didn't show up in studlistACTIVE or studlistMARGINAL, and therefore completed no learning
        %%%%tasks 
        AchievementData.Q(ii) = 0 ;
        AchievementData.T(ii) = 0 ;
        g2.minspent(g2.student_id==AchievementData.id(ii)) = AchievementData.T(ii) ;  %g2 contains a different and less relevant measure of time spent, while ProductionDataNonCum is the more relevant and reliable measure, based on accumulating task-by-task pageload times.
    else %%%%If we make it here, the student completed at least 1 learning task
        AchievementData.Q(ii)      = g2.total_pass(g2.student_id==AchievementData.id(ii)) ; 
        AchievementData.T(ii)      = sum(ProductionDataNonCum.t(ProductionDataNonCum.sid==AchievementData.id(ii))) + unfinished.time(unfinished.sid==AchievementData.id(ii)) ; 
        g2.minspent(g2.student_id==AchievementData.id(ii)) = AchievementData.T(ii) ;  %g2 contains a different and less relevant measure of time spent, while ProductionDataNonCum is the more relevant and reliable measure, based on accumulating task-by-task pageload times.
        AchievementData.Tunfinished(ii) = unfinished.time(unfinished.sid==AchievementData.id(ii)) ; 
        %%%%The following will contain work time on the first unit of success (which is NOT the same
        %%%%as total time for marginal students; this also included "unfinished" time...).  This is
        %%%%the same measure as EPointEst.FirstTime
        AchievementData.Tfirstunit(ii) = EPointEst.FirstTime(EPointEst.studlistFULL==AchievementData.id(ii)) ;
        EPointEstindx              = find(studlistFULL==AchievementData.id(ii)) ;  %%%%This is the index within EPointEst.QE/T/ThetaE, which contains marginal and active students
        AchievementData.ThetaE(ii) = EPointEst.ThetaE(EPointEstindx) ; 
    end
    
    %%%%In this loop we fill in the base and marginal wage numbers:
    if g2.contract1(g2.student_id==AchievementData.id(ii))==1
        AchievementData.contract(ii) = 1 ;
        AchievementData.b(ii) = 15 ;
        AchievementData.w(ii) = 0.75 ;
    elseif g2.contract2(g2.student_id==AchievementData.id(ii))==1
        AchievementData.contract(ii) = 2 ;
        AchievementData.b(ii) = 10 ;
        AchievementData.w(ii) = 1 ;
    else
        AchievementData.contract(ii) = 3 ;
        AchievementData.b(ii) = 5 ;
        AchievementData.w(ii) = 1.25 ;
    end
end

%%%%This last step is a small housecleaning measure: For anyone in the marginal list who's mean time
%%%%to completion is lower than a standard deviation below the active mean, we assume their first
%%%%(and only) success was tainted too much by luck, and we drop them from the marginal group and
%%%%reset the values of Q, T, Tfirstunit, and ThetaE.
marginaldropcutoff = mean(AchievementData.T(AchievementData.Q>=2)./AchievementData.Q(AchievementData.Q>=2))...
                          - std(AchievementData.T(AchievementData.Q>=2)./AchievementData.Q(AchievementData.Q>=2)) ;
studlistMARGINAL_droplist = [] ;
for ii=1:length(studlistMARGINAL)
    if AchievementData.T(AchievementData.id==studlistMARGINAL(ii))<marginaldropcutoff
        disp('DROPPING ONE NOW!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
        AchievementData.Q(AchievementData.id==studlistMARGINAL(ii)) = 0 ;
        AchievementData.T(AchievementData.id==studlistMARGINAL(ii)) = 0 ;
        AchievementData.Tunfinished(AchievementData.id==studlistMARGINAL(ii)) = 0 ;
        AchievementData.Tfirstunit(AchievementData.id==studlistMARGINAL(ii)) = 0 ;
        AchievementData.ThetaE(AchievementData.id==studlistMARGINAL(ii)) = nan ;
        
        studlistMARGINAL_droplist = [studlistMARGINAL_droplist; studlistMARGINAL(ii)] ;  %#ok  %%%%This will be used later to expunge the corresponding student ids from studlistMARGINAL and studlistFULL
        
        EPointEst.ThetaEMarginal(EPointEst.marginallist==studlistMARGINAL(ii)) = [] ;
        EPointEst.marginallist(EPointEst.marginallist==studlistMARGINAL(ii)) = [] ;
        EPointEst.ThetaE(EPointEst.studlistFULL==studlistMARGINAL(ii)) = [] ;
        EPointEst.FirstTime(EPointEst.studlistFULL==studlistMARGINAL(ii)) = [] ;
        EPointEst.studlistFULL(EPointEst.studlistFULL==studlistMARGINAL(ii)) = [] ;

        g2.total_attempt(g2.student_id==studlistMARGINAL(ii)) = 0 ;
        g2.total_pass(g2.student_id==studlistMARGINAL(ii)) = 0 ;
        g2.timing_pagesubmit(g2.student_id==studlistMARGINAL(ii)) = 0 ;
        g2.minspent(g2.student_id==studlistMARGINAL(ii)) = 0 ;
        g2.succrat(g2.student_id==studlistMARGINAL(ii)) = nan ;
    end
end
for ii=1:length(studlistMARGINAL_droplist)
    studlistMARGINAL(studlistMARGINAL==studlistMARGINAL_droplist(ii)) = [] ;
    studlistFULL(studlistFULL==studlistMARGINAL_droplist(ii)) = [] ;
end


%%%%DEMOGRAPHIC VARIABLES:
AchievementData.fem                 = (g2.gender=="f") ; %%%%This is a dummy for whether the kid is female
AchievementData.fem(g2.gender==".") = -1  ; %%%%if the observation is missing, we set its value to -1
AchievementData.white               = g2.white ;
AchievementData.asian               = g2.asian ;  AchievementData.asian(isnan(AchievementData.asian)) = 0 ;
AchievementData.black               = g2.black ;
AchievementData.hisp                = g2.hispanic ;
%%%%The following variables are for the census block group (cbg) where the student resides:
AchievementData.mean_inc_cbg             = g2.mean_inc_cbg ; %%%%This is mean household income
 
%%%%MATH EXAM VARIABLES:
AchievementData.pretest         = g2.pre_rawscore ;  %%%%This is the student's pre-test score (out of 36)
AchievementData.preblank        = g2.totalpreblank ; %%%%This is how many questions on te pre-test were left blank
AchievementData.posttest        = g2.post_rawscore ;  %%%%This is the student's post-test score (out of 36)
AchievementData.postblank       = g2.totalpostblank ;  %%%%This is how many questions on te post-test were left blank
AchievementData.grade5          = g2.grade5 ;  %%%%This is a dummy for the student's grade-level cohort  Group 2 has only 5th and 6th graders.

%%%%DEVICE AVAILABILITY VARIABLES:
AchievementData.mobile         = g2.mobile ;  %%%%This is a count of how many times the student connected to ChicagoMathGame.org from a smartphone
AchievementData.computer       = g2.computer ;  %%%%This is a count of how many times the student connected to ChicagoMathGame.org from a desktop or laptop computer
AchievementData.tablet         = g2.tablet ;  %%%%This is a count of how many times the student connected to ChicagoMathGame.org from a tablet
AchievementData.ngamingsys     = g2.n_gaming ;  %%%%This is a count of how many video gaming systems the kid has at home (being Xbox, Playstation, Nintendo, or Other)
AchievementData.gaming_weekday = g2.gaming_weekday ;  %%%%This is a dummy variable for whether the kid's parents allow him/her to play video games on weekdays.
AchievementData.no_internet    = g2.no_internet ;  %%%%This is a dummy for whether the kid has an internet connection at home.
AchievementData.weekendnet     = g2.weekend_internet ; %%%%A value of 1 indicates that the kid has internet at home, but is not allowed to use it for recreational purposes on weekdays
AchievementData.weekdaynet1    = g2.onehour_internet_weekday ; %%%%A value of 1 indicates that the kid has internet at home, and is allowed to use it for recreational purposes for one hour daily on weekdays
AchievementData.weekdaynet2    = g2.twohour_internet_weekday ; %%%%A value of 1 indicates that the kid has internet at home, and is allowed to use it for recreational purposes for two hours daily on weekdays
AchievementData.weekdaynet3    = g2.twoplus_internet_weekday ; %%%%A value of 1 indicates that the kid has internet at home, and is allowed to use it for recreational purposes for more than two hours daily on weekdays


%%%%SELF-REPORTED HELPER VARIABLES:
AchievementData.helpmom     = (g2.helpcmg_mom>0 | g2.helphw_mom>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their mother helps them with academic work when needed.
AchievementData.helpdad     = (g2.helpcmg_dad>0 | g2.helphw_dad>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their father helps them with academic work when needed.
AchievementData.helpsis     = (g2.helpcmg_sis>0 | g2.helphw_sis>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their sister helps them with academic work when needed.
AchievementData.helpbro     = (g2.helpcmg_bro>0 | g2.helphw_bro>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their brother helps them with academic work when needed.
AchievementData.helpgrand   = (g2.helpcmgoth_grand>0 | g2.helphwoth_grand>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their grandparent helps them with academic work when needed.
AchievementData.helptut     = (g2.helpcmgoth_tutor>0 | g2.helphwoth_tutor>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their tutor helps them with academic work.
AchievementData.helppal     = (g2.helpcmgoth_friend>0 | g2.helphwoth_friend>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that their friend helps them with academic work when needed.
AchievementData.helpoth     = (g2.helpcmgoth_another>0 | g2.helphwoth_another>0) ; %%%%This is a binary for whether the student reported on either the pre-survey or post-survey that someone else not mentioned above helps them with academic work when needed.
AchievementData.helpadult   = g2.helpadult ;  %%%%This is a count of how many adult helpers the student self-reported
AchievementData.helppeer    = g2.helppeer;  %%%%This is a count of how many peer helpers the student self-reported

%%%%TIME USE VARIABLES:
AchievementData.time_social_unstr = g2.time_social_unstructured ;  %%%%This is the self-reported amount of time in a typical day (averaging over weekday and weekend responses) that the student spends in NON-adult-supervised time with friends
AchievementData.time_social_str   = g2.time_social_structured ;  %%%%This is the self-reported amount of time in a typical day (averaging over weekday and weekend responses) that the student spends in adult-supervised time with friends
AchievementData.time_screen       = g2.time_screen ;  %%%%This is the self-reported amount of time on a typical day (averaging over weekday and weekend responses) that the student spends in front of an electronic screen for recreational purposes
AchievementData.time_study        = g2.time_study ;  %%%%This is the self-reported amount of time on a typical day (averaging over weekday and weekend responses) that the student spends on academic study
AchievementData.freq_miss_hmwk    = g2.freq_miss_hmwk ;  %%%%This is the self-reported fraction of the time that a student fails to turn in math homework.  "NEVER"=0, "SOMETIMES"=0.25, "ABOUT HALF OF THE TIME"=0.5, "MOST OF THE TIME"=0.75, "ALMOST ALWAYS"=1
AchievementData.sports            = g2.sports_d ; %%%%This is a dummy for whether the kid is involved in any formal sports activities
AchievementData.music             = g2.music_d ; %%%%This is a dummy for whether the kid is involved in any formal music lessons/activities
AchievementData.otheractivity     = g2.other_d ; %%%%This is a dummy for whether the kid is involved in any other formal clubs/activities outside of sports or music

%%%%ACADEMIC ATTITUDE VARIABLES: 
AchievementData.math_favorite     = g2.math_favorite ;  %%%%This is a dummy for whether the student specified MATH as their favorite out of the following list: math, science, english, history/social studies, foreign language studies
AchievementData.science_favorite  = g2.science_favorite ;  %%%%This is a dummy for whether the student specified SCIENCE as their favorite out of the following list: math, science, english, history/social studies, foreign language studies
AchievementData.math_least        = g2.math_least ;  %%%%This is a dummy for whether the student specified MATH as their LEAST favorite out of the following list: math, science, english, history/social studies, foreign language studies
AchievementData.science_least     = g2.science_least ;  %%%%This is a dummy for whether the student specified SCIENCE as their LEAST favorite out of the following list: math, science, english, history/social studies, foreign language studies
AchievementData.STEMfav           = g2.math_favorite + g2.science_favorite ; %%%%This is a dummy for whether the student self reported EITHER math OR science as their favorite subject from the list above.
%%%%The following are two measures of how extrinsically motivated the student is.  There were 4 
%%%%survey questions (2 pre and 2 post) asking what their biggest motivation was for doing math work 
%%%%and what their second biggest motivation was.  Choices included intrinsic and extrinsic 
%%%%responses, and a "none of the above" option.  This measure assigns 1 point for each biggest 
%%%%motivation response and 0.5 point for each second biggest motivation response.                                                                
AchievementData.extrinsic         = g2.mot_extr ;                                                            
AchievementData.intrinsic         = g2.mot_intr ; 

%%%%SCHOOL DISTRICT FIXED EFFECTS:
AchievementData.Dist1 = (g2.district_id==1) ; %%%%This is a dummy for whether the kid attends District 1
AchievementData.Dist2 = (g2.district_id==3) ; %%%%This is a dummy for whether the kid attends District 2
AchievementData.Dist3 = (g2.district_id==4) ; %%%%This is a dummy for whether the kid attends District 3

                          


format shortg

run AAA_FILE02_DescriptiveAnalyses.m
%%

disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGE 2 ESTIMATION:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')

run AAA_FILE03_ThetaLDriver.m  








AchievementData.min    = zeros(size(AchievementData.ThetaE)) ;
AchievementData.ThetaL = nan(size(AchievementData.ThetaE)) ;
for ii=1:length(AchievementData.ThetaE)
    if AchievementData.black(ii)==1 | AchievementData.hisp(ii)==1 %#ok
        AchievementData.min(ii) = 1 ;
    end
    if AchievementData.Q(ii)>=2 && AchievementData.Q(ii)<80
        tempindx = find(LPointEst.student_id==AchievementData.id(ii)) ;
        AchievementData.ThetaL(ii) = LPointEst.ThetaLACT(tempindx) ;
    elseif AchievementData.Q(ii)==80 
        tempcon = AchievementData.contract(ii) ;
        Qwts    = LPointEst.Q80freq(:,2,tempcon) ;
        Qwts    = Qwts(~isnan(Qwts)) ;
        Qwts    = Qwts/sum(Qwts) ;
        AchievementData.ThetaL(ii) = sum( Qwts.*LPointEst.ThetaL80(LPointEst.ThetaL80(:,1)==AchievementData.id(ii),2:1+sum(~isnan(LPointEst.Q80freq(:,2,tempcon))))' ) ;
    end
end


%%
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGES 1 & 2 BOOTSTRAP ESTIMATION:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%       NOW RUNNING THE BOOTSTRAP:       %%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% BOOTSTRAP CODE: THIS SEGMENT OF CODE RUNS THE RE-SAMPLING SCHEME.  THE NEXT FILE CALLED BELOW, 
%%%% AAA_FILE04_BootstrapDriver_MultiMethod7.m ACTUALLY COMPUTES BOOTSTRAP ESTIMATES FOR EACH
%%%% BOOTSTRAP SAMPLE.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%When we randomized we balanced on race, gender, grade and school district, and we also 
%%%%stratified on pre-test grades).  Thus, here we do a block bootstrap routine where we first bin
%%%%the data by treatment-race-gender-grade-school district bins, and then we sample with 
%%%%replacement from those bins.
clear bin_*


indx                  = find( (AchievementData.black==1|AchievementData.hisp==1) ...
                               & AchievementData.fem==1 & AchievementData.Dist2==1 ...
                               & AchievementData.grade5==1 ) ; 
bin_fem_min_Dist2_5 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( (AchievementData.black==1|AchievementData.hisp==1) ...
                               & AchievementData.fem==1 & AchievementData.Dist2==1 ...
                               & AchievementData.grade5==0 ) ; 
bin_fem_min_Dist2_6 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( ~(AchievementData.black==1|AchievementData.hisp==1) ...
                                & AchievementData.fem==1 & AchievementData.Dist2==1 ...
                                & AchievementData.grade5==1 ) ; 
bin_fem_non_Dist2_5 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( ~(AchievementData.black==1|AchievementData.hisp==1) ...
                                & AchievementData.fem==1 & AchievementData.Dist2==1 ...
                                & AchievementData.grade5==0 ) ; 
bin_fem_non_Dist2_6 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( AchievementData.fem==1 & AchievementData.Dist3==1 ...
                              & AchievementData.grade5==1 ) ; 
bin_fem_Dist3_5      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( AchievementData.fem==1 & AchievementData.Dist3==1 ...
                              & AchievementData.grade5==0 ) ; 
bin_fem_Dist3_6      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( AchievementData.fem==1 & AchievementData.Dist1==1 ...
                              & AchievementData.grade5==1 ) ; 
bin_fem_Dist1_5      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                  = find( AchievementData.fem==1 & AchievementData.Dist1==1 ...
                              & AchievementData.grade5==0 ) ; 
bin_fem_Dist1_6      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ;  

indx                   = find( (AchievementData.black==1|AchievementData.hisp==1) ...
                                & AchievementData.fem==0 & AchievementData.Dist2==1 ...
                                & AchievementData.grade5==1 ) ; 
bin_male_min_Dist2_5 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( (AchievementData.black==1|AchievementData.hisp==1) ...
                                & AchievementData.fem==0 & AchievementData.Dist2==1 ...
                                & AchievementData.grade5==0 ) ; 
bin_male_min_Dist2_6 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( ~(AchievementData.black==1|AchievementData.hisp==1) ...
                                 & AchievementData.fem==0 & AchievementData.Dist2==1 ...
                                 & AchievementData.grade5==1 ) ; 
bin_male_non_Dist2_5 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( ~(AchievementData.black==1|AchievementData.hisp==1) ...
                                 & AchievementData.fem==0 & AchievementData.Dist2==1 ...
                                 & AchievementData.grade5==0 ) ; 
bin_male_non_Dist2_6 = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( AchievementData.fem==0 & AchievementData.Dist3==1 ...
                               & AchievementData.grade5==1 ) ; 
bin_male_Dist3_5      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( AchievementData.fem==0 & AchievementData.Dist3==1 ...
                               & AchievementData.grade5==0 ) ; 
bin_male_Dist3_6      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( AchievementData.fem==0 & AchievementData.Dist1==1 ...
                               & AchievementData.grade5==1 ) ; 
bin_male_Dist1_5      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 

indx                   = find( AchievementData.fem==0 & AchievementData.Dist1==1 ...
                               & AchievementData.grade5==0 ) ; 
bin_male_Dist1_6      = [AchievementData.id(indx) AchievementData.pretest(indx) AchievementData.contract(indx) AchievementData.Q(indx)] ; 


rng(round(exp(1)*1000000)); 
BootstrapSims = 6400 ;     
BSSamp = cell(BootstrapSims,1) ;
binnames = {'bin_fem_Dist3_5';'bin_fem_Dist3_6';'bin_fem_Dist1_5';'bin_fem_Dist1_6';'bin_fem_min_Dist2_5';'bin_fem_min_Dist2_6';'bin_fem_non_Dist2_5';'bin_fem_non_Dist2_6';'bin_male_Dist3_5';'bin_male_Dist3_6';'bin_male_Dist1_5';'bin_male_Dist1_6';'bin_male_min_Dist2_5';'bin_male_min_Dist2_6';'bin_male_non_Dist2_5';'bin_male_non_Dist2_6'} ;
for s=1:BootstrapSims
    BSSamp{s,1} = nan(length(AchievementData.id),1) ;
    for kk=1:length(binnames)
        cmd = strcat('temp=',binnames{kk},';') ;
        eval(cmd);
        temp1 = temp(temp(:,3)==1,:) ; 
        temp2 = temp(temp(:,3)==2,:) ; 
        temp3 = temp(temp(:,3)==3,:) ; 
        
        N_1       = length(temp1(:,1)) ;
        N_2       = length(temp2(:,1)) ;
        N_3       = length(temp3(:,1)) ;
        N         = max([N_1;N_2;N_3]) ;
        tempsamp1 = nan(N,1) ;
        tempsamp2 = nan(N,1) ;
        tempsamp3 = nan(N,1) ;
        sampindx  = ceil(length(temp(:,1))*rand(N,1)) ;
        for jj=1:N
            if temp(sampindx(jj),3)==1 %%%%That means this person was in contract 1
                tempsamp1(jj) = temp(sampindx(jj),1) ;
                tempscore     = temp(sampindx(jj),2) ;  %%%%This is the person's pre-test score.
                %%%%Below we select two people from the other two contract groups in the same bin.
                %%%%In order to mimic the balancing on pre-test scores, those two people are chosen
                %%%%to have scores as close as possible to the first person's score.  Below we do a
                %%%%robust form of nearest-neighbor interpolation that randomly breaks ties if they
                %%%%arise.
                temptemp2     = abs(temp2(:,2)-tempscore) ;
                temptemp2     = temp2(temptemp2==min(temptemp2),1) ;
                if length(temptemp2)==1
                    tempsamp2(jj) = temptemp2 ;
                else
                    tempsamp2(jj) = temptemp2( ceil(rand(1,1)*length(temptemp2)) ) ;
                end
                temptemp3     = abs(temp3(:,2)-tempscore) ;
                temptemp3     = temp3(temptemp3==min(temptemp3),1) ;
                if length(temptemp3)==1
                    tempsamp3(jj) = temptemp3 ;
                else
                    tempsamp3(jj) = temptemp3( ceil(rand(1,1)*length(temptemp3)) ) ;
                end
            elseif temp(sampindx(jj),3)==2 %%%%That means this person was in contract 2
                tempsamp2(jj) = temp(sampindx(jj),1) ;
                tempscore     = temp(sampindx(jj),2) ;  %%%%This is the person's pre-test score.
                %%%%Below we select two people from the other two contract groups in the same bin.
                %%%%In order to mimic the balancing on pre-test scores, those two people are chosen
                %%%%to have scores as close as possible to the first person's score.  Below we do a
                %%%%robust form of nearest-neighbor interpolation that randomly breaks ties if they
                %%%%arise.
                temptemp1     = abs(temp1(:,2)-tempscore) ;
                temptemp1     = temp1(temptemp1==min(temptemp1),1) ;
                if length(temptemp1)==1
                    tempsamp1(jj) = temptemp1 ;
                else
                    tempsamp1(jj) = temptemp1( ceil(rand(1,1)*length(temptemp1)) ) ;
                end
                temptemp3     = abs(temp3(:,2)-tempscore) ;
                temptemp3     = temp3(temptemp3==min(temptemp3),1) ;
                if length(temptemp3)==1
                    tempsamp3(jj) = temptemp3 ;
                else
                    tempsamp3(jj) = temptemp3( ceil(rand(1,1)*length(temptemp3)) ) ;
                end
            elseif temp(sampindx(jj),3)==3 %%%%That means this person was in contract 3
                tempsamp3(jj) = temp(sampindx(jj),1) ;
                tempscore     = temp(sampindx(jj),2) ;  %%%%This is the person's pre-test score.
                %%%%Below we select two people from the other two contract groups in the same bin.
                %%%%In order to mimic the balancing on pre-test scores, those two people are chosen
                %%%%to have scores as close as possible to the first person's score.  Below we do a
                %%%%robust form of nearest-neighbor interpolation that randomly breaks ties if they
                %%%%arise.
                temptemp1     = abs(temp1(:,2)-tempscore) ;
                temptemp1     = temp1(temptemp1==min(temptemp1),1) ;
                if length(temptemp1)==1
                    tempsamp1(jj) = temptemp1 ;
                else
                    tempsamp1(jj) = temptemp1( ceil(rand(1,1)*length(temptemp1)) ) ;
                end
                temptemp2     = abs(temp2(:,2)-tempscore) ;
                temptemp2     = temp2(temptemp2==min(temptemp2),1) ;
                if length(temptemp2)==1
                    tempsamp2(jj) = temptemp2 ;
                else
                    tempsamp2(jj) = temptemp2( ceil(rand(1,1)*length(temptemp2)) ) ;
                end
            end
        end
        last                         = find(isnan(BSSamp{s,1}(:)),1,'first') - 1 ;
        BSSamp{s,1}(last+1:last+N_1) = tempsamp1(1:N_1) ;
        last                         = last + N_1 ;
        BSSamp{s,1}(last+1:last+N_2) = tempsamp2(1:N_2) ;
        last                         = last + N_2 ;
        BSSamp{s,1}(last+1:last+N_3) = tempsamp3(1:N_3) ;
    end
    BSSamp{s,1} = sort(BSSamp{s,1}) ;
end

save('PointEstimates_MSM7.mat','EPointEst','LPointEst','AchievementData','g2','ProductionDataNonCum','BSSamp','BootstrapSims')




%%

disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGES 1 & 2 PRE-COMPUTED BOOTSTRAP ESTIMATES:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  NOTES ON COMPUTATION TME: THE BOOTSTRAP HERE IS QUITE TIME CONSUMING FOR TWO REASONS:
%%%%%%%  OBJECTIVE FUNCTION EVALUATIONS TAKE A LOT OF TIME, AND THE ESTIMATOR MUST BE RUN MANY TIMES
%%%%%%%  OVER.  THAT BEING THE CASE, THE OBJECTIVE FUNCTION EVALUATIONS ARE PARALLELIZED.  ALTHOUGH
%%%%%%%  THIS MEANS THAT THE BOOTSTRAPS CANNOT ALSO BE DIRECTLY PARALLELIZED, THIS CAN BE DONE IN
%%%%%%%  BATCHES.  THE OPTIMAL NUMBER OF WORKERS TO DEVOTE TO THE PARALLEL OBJECTIVE FUNCTION
%%%%%%%  COMPONENT IS AROUND 20-30 (GIVEN SOLVER METHODS USED), BUT MULTIPLE INSTANCES OF MATLAB
%%%%%%%  CAN BE OPEN AT THE SAME TIME ON A SINGLE MULTI-PROCESSOR COMPUTER, WITH EACH INSTANCE
%%%%%%%  HAVING THE NUMBER OF WORKERS ABOVE, UNTIL ALL CORES ON THE COMPUTER ARE MAXED OUT.  THIS
%%%%%%%  ALLOWS FOR PARALLELIZATION BOTH IN THE OBJECTIVE FUNCTION EVALS AND OF THE BOOTSTRAPS TOO.
%%%%%%%  THIS IS HOW BOOTSTRAPS FOR THE PAPER WERE PERFORMED, AND THIS IS WHY BootstrapSims IS SO
%%%%%%%  LARGE, EVEN THOUGH THE TARGET IS ONLY 400 BOOTSTRAPS.  TO SEPARATE THE BOOTSTRAPS INTO
%%%%%%%  BATCHES, SIMPLY ALTER THE "start" AND "stop" VALUES BELOW WITHIN EACH MATLAB INSTANCE.
%%%%%%%  E.G., ONE INSTANCE COULD USE "start=1; stop=50" WHILE A SECOND COULD USE  "start=51;
%%%%%%%  stop=100", AND SO ON (MAKING SURE THAT stop<=6400)...  THE FUNCTION PROCESSBOOTSTRAPS.M 
%%%%%%%  CAN COMBINE AND PROCESS THE BOOTSTRAP BATCHES, POST-ESTIMATION. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all %#ok
RunBootstrap = 0 ;
if RunBootstrap==1
    start  = 1 %#ok
    finish = 100 %#ok
    run AAA_FILE04_BootstrapDriver_MultiMethod7
else
    run ProcessBootstraps
end


      


basepay        = 15*(AchievementData.contract==1) + 10*(AchievementData.contract==2) + 5*(AchievementData.contract==3) ;
margpay        = 0.75*(AchievementData.contract==1) + 1*(AchievementData.contract==2) + 1.25*(AchievementData.contract==3) ;
wageratehr     = (AchievementData.Q(AchievementData.Q>=2).*margpay(AchievementData.Q>=2)+basepay(AchievementData.Q>=2))./(AchievementData.T(AchievementData.Q>=2)/60)  ;  
margwageratehr = (AchievementData.Q(AchievementData.Q>=2).*margpay(AchievementData.Q>=2))./(AchievementData.T(AchievementData.Q>=2)/60)  ;   





wtsE              = 1./(AchievementData.lThetaEse.^2) ; 
wtsE(isnan(wtsE)) = min(wtsE) ; 
wtsE              = length(wtsE)*wtsE/sum(wtsE,'omitnan') ;
[mean(wtsE,'omitnan') std(wtsE,'omitnan') min(wtsE) max(wtsE)] %#ok
wtsL              = 1./(AchievementData.lThetaLse.^2) ; 
wtsL(isnan(wtsL)) = min(wtsL) ; 
wtsL              = (length(wtsL)/sum(wtsL,'omitnan'))*wtsL ;
wtsL(wtsL<min(wtsE)) = min(wtsE) ;
wtsL(wtsL>max(wtsE)) = max(wtsE) ;
wtsL                 = (length(wtsL)/sum(wtsL,'omitnan'))*wtsL ;
[mean(wtsL,'omitnan') std(wtsL,'omitnan') min(wtsL) max(wtsL)] %#ok
[sum(isnan(wtsE)) sum(isnan(wtsL))] %#ok
wts = wtsE ;  %%%%Should just use the weights on ThetaE, since those are directly determined by
%%%%panel length, and they indirectly determine std errors on TheteL




domE = linspace(EPointEst.cutoffs(1),EPointEst.cutoffs(end),100)' ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%NOW COMPUTING THE SELECTION THRESHOLD FOR CONTRACT GROUP 1:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x      = log(AchievementData.ThetaE(AchievementData.contract==1&AchievementData.Q>=2));
y      = log(AchievementData.ThetaL(AchievementData.contract==1&AchievementData.Q>=2));
k      = convhull(x,y);
k(end) = [] ;
%%%%Next we convert the convex hull into a convex, non-increasing upper-bound by selecting all
%%%%members of the convex hull boundary (x(k),y(k)) that have no other neighbors in the set
%%%%(x,y) to their northwest.
upperboundflag = zeros(length(k),1);
for ii=1:length(k)
    if sum( (x>=x(k(ii)) & y>y(k(ii))) | (x>x(k(ii)) & y>=y(k(ii))) )==0
        upperboundflag(ii)=1;
    end
end
threshold1 = sortrows([x(k(upperboundflag==1)) y(k(upperboundflag==1))]);
threshold1 = exp( interp1(threshold1(:,1),threshold1(:,2),log(domE),'linear','extrap') ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%NOW COMPUTING THE SELECTION THRESHOLD FOR CONTRACT GROUP 2:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x      = log(AchievementData.ThetaE(AchievementData.contract==2&AchievementData.Q>=2));
y      = log(AchievementData.ThetaL(AchievementData.contract==2&AchievementData.Q>=2));
k      = convhull(x,y);
k(end) = [] ;
%%%%Next we convert the convex hull into a convex, non-increasing upper-bound by selecting all
%%%%members of the convex hull boundary (x(k),y(k)) that have no other neighbors in the set
%%%%(x,y) to their northwest.
upperboundflag = zeros(length(k),1);
for ii=1:length(k)
    if sum( (x>=x(k(ii)) & y>y(k(ii))) | (x>x(k(ii)) & y>=y(k(ii))) )==0
        upperboundflag(ii)=1;
    end
end
threshold2 = sortrows([x(k(upperboundflag==1)) y(k(upperboundflag==1))]);
threshold2 = exp( interp1(threshold2(:,1),threshold2(:,2),log(domE),'linear','extrap') ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%NOW COMPUTING THE SELECTION THRESHOLD FOR CONTRACT GROUP 3:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x      = log(AchievementData.ThetaE(AchievementData.contract==3&AchievementData.Q>=2));
y      = log(AchievementData.ThetaL(AchievementData.contract==3&AchievementData.Q>=2));
k      = convhull(x,y);
k(end) = [] ;
%%%%Next we convert the convex hull into a convex, non-increasing upper-bound by selecting all
%%%%members of the convex hull boundary (x(k),y(k)) that have no other neighbors in the set
%%%%(x,y) to their northwest.
upperboundflag = zeros(length(k),1);
for ii=1:length(k)
    if sum( (x>=x(k(ii)) & y>y(k(ii))) | (x>x(k(ii)) & y>=y(k(ii))) )==0
        upperboundflag(ii)=1;
    end
end
threshold3 = sortrows([x(k(upperboundflag==1)) y(k(upperboundflag==1))]);
threshold3 = exp( interp1(threshold3(:,1),threshold3(:,2),log(domE),'linear','extrap') ) ;
LPointEst.domE = domE ;
LPointEst.threshold1 = threshold1 ;
LPointEst.threshold2 = threshold2 ;
LPointEst.threshold3 = threshold3 ;




AchievementData.weekdaynet        = 0 + 1*AchievementData.weekdaynet1 + 2*AchievementData.weekdaynet2 + 3*AchievementData.weekdaynet3 ;
AchievementData.mathpref          = -1*AchievementData.math_least + 1*AchievementData.math_favorite ;
AchievementData.mobilefrac        = AchievementData.mobile./(AchievementData.mobile+AchievementData.computer+AchievementData.tablet) ; 
AchievementData.mobilefrac(AchievementData.Q<2) = nan ;
AchievementData.computerfrac      = AchievementData.computer./(AchievementData.mobile+AchievementData.computer+AchievementData.tablet) ;
AchievementData.computerfrac(AchievementData.Q<2) = nan ; 
AchievementData.tabletfrac        = AchievementData.tablet./(AchievementData.mobile+AchievementData.computer+AchievementData.tablet) ; 
AchievementData.tabletfrac(AchievementData.Q<2) = nan ;
AchievementData.time_screen_pre   = g2.time_screen_pre ;  %%%%This is the self-reported amount of time on a typical day (averaging over weekday and weekend responses) that the student spends in front of an electronic screen for recreational purposes
AchievementData.understand_better = g2.understand_better ; AchievementData.understand_better(isnan(g2.understand_better)) = nan ;
AchievementData.understand_lbetter= g2.understand_littlebetter ; AchievementData.understand_better(isnan(g2.understand_littlebetter)) = nan ;
AchievementData.enjoy_better      = g2.enjoy_better ;  
AchievementData.enjoy_lbetter     = g2.enjoy_littlebetter ;  
AchievementData.miss_hmwk_some    = 1*(g2.freq_miss_hmwk==0.25) ; AchievementData.miss_hmwk_some(isnan(g2.freq_miss_hmwk)) = nan ;
AchievementData.miss_hmwk_usually = 1*(g2.freq_miss_hmwk>=0.5) ; AchievementData.miss_hmwk_usually(isnan(g2.freq_miss_hmwk)) = nan ;
AchievementData.bigfam              = g2.bigfam ;
AchievementData.involved_parent     = g2.involved_parent ; 
AchievementData.educator_parent     = 1*(g2.edu_occupation_reson==1 | g2.edu_occupcation_partner==1) ;
AchievementData.educator_parent(isnan(g2.edu_occupation_reson) & isnan(g2.edu_occupcation_partner)) = nan ; 
AchievementData.STEM_occ_parent     = 1*(g2.STEM_occupation_respon==1 | g2.STEM_occupation_partner==1) ;
AchievementData.STEM_occ_parent(isnan(g2.STEM_occupation_respon) & isnan(g2.STEM_occupation_partner)) = nan  ;
AchievementData.qualified_parent    = 1*(AchievementData.educator_parent==1 | AchievementData.STEM_occ_parent==1) ; AchievementData.qualified_parent(isnan(AchievementData.educator_parent)&isnan(AchievementData.STEM_occ_parent)) = nan ;
AchievementData.disposable_inc      = g2.disposable_inc ; 
AchievementData.oldest_child        = g2.oldest_child ; AchievementData.oldest_child(g2.only_child==1)=1 ;  %%%NOTE: "oldest" means either youngest or an only child.
AchievementData.middle_child        = g2.middle_child ; 
AchievementData.youngest_child      = g2.youngest_child; 
AchievementData.oldest_child(g2.oldest_child==0&g2.middle_child==0&g2.youngest_child==0&g2.only_child==0) = nan ;
AchievementData.middle_child(g2.oldest_child==0&g2.middle_child==0&g2.youngest_child==0&g2.only_child==0) = nan ;
AchievementData.youngest_child(g2.oldest_child==0&g2.middle_child==0&g2.youngest_child==0&g2.only_child==0) = nan ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%                 FIGURE 5:                 %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure; 
subplot(1,2,1); 
hold on; 
ecdf(AchievementData.Q(AchievementData.black==1|AchievementData.hisp==1)); 
ecdf(AchievementData.Q(AchievementData.white==1|AchievementData.asian==1)); 
legend('BLACK/HISPANIC','WHITE/ASIAN','Location','southeast','FontWeight','bold','FontSize',28);
box on; 
grid on; 
xlabel('WEBSITE LEARNING TASK COMPLETION','FontWeight','bold','FontSize',28); 
ylabel('CDFs','FontWeight','bold','FontSize',28); 
ylim([0.4,1]); 
subplot(1,2,2); 
hold on; 
ecdf(AchievementData.Q(AchievementData.fem==0)); 
ecdf(AchievementData.Q(AchievementData.fem==1)); 
legend('MALE','FEMALE','Location','southeast','FontWeight','bold','FontSize',28); 
box on; 
grid on; 
xlabel('WEBSITE LEARNING TASK COMPLETION','FontWeight','bold','FontSize',28); 
ylabel('CDFs','FontWeight','bold','FontSize',28); 
ylim([0.4,1])

%%
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGE 3 ESTIMATION:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%         TABLE 7: COVARIATE DESCRIPTIVE STATS')
disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%--------------------------------------------------------------------------------------')

MuSigFMO_Fem = [mean(AchievementData.fem,'omitnan') std(AchievementData.fem,'omitnan') mean(isnan(AchievementData.fem))] %#ok
MuSigFMO_Male = [mean(1-AchievementData.fem,'omitnan') std(1-AchievementData.fem,'omitnan') mean(isnan(1-AchievementData.fem))] %#ok
MuSigFMO_Blk = [mean(AchievementData.black,'omitnan') std(AchievementData.black,'omitnan') mean(isnan(AchievementData.black))] %#ok
MuSigFMO_Hsp = [mean(AchievementData.hisp,'omitnan') std(AchievementData.hisp,'omitnan') mean(isnan(AchievementData.hisp))] %#ok
MuSigFMO_WhtAsian = [mean(1-AchievementData.black-AchievementData.hisp,'omitnan') std(1-AchievementData.black-AchievementData.hisp,'omitnan') mean(isnan(1-AchievementData.black-AchievementData.hisp))] %#ok
MuSigFMO_Grade5 = [mean(AchievementData.grade5,'omitnan') std(AchievementData.grade5,'omitnan') mean(isnan(AchievementData.grade5))] %#ok
MuSigFMO_NbhdInc = [mean(AchievementData.mean_inc_cbg/1000,'omitnan') std(AchievementData.mean_inc_cbg/1000,'omitnan') mean(isnan(AchievementData.mean_inc_cbg))] %#ok
MuSigFMO_NbhdIncBlk = [mean(AchievementData.mean_inc_cbg(AchievementData.black==1)/1000,'omitnan')] %#ok
MuSigFMO_NbhdIncHsp = [mean(AchievementData.mean_inc_cbg(AchievementData.hisp==1)/1000,'omitnan')] %#ok
MuSigFMO_NbhdIncWA = [mean(AchievementData.mean_inc_cbg(AchievementData.black==0&AchievementData.hisp==0)/1000,'omitnan')] %#ok
MuSigFMO_NbhdIns = [mean(AchievementData.minor_no_healthins_cbg,'omitnan') std(AchievementData.minor_no_healthins_cbg,'omitnan') mean(isnan(AchievementData.minor_no_healthins_cbg))] %#ok
MuSigFMO_NbhdInsBlk = [mean(AchievementData.minor_no_healthins_cbg(AchievementData.black==1),'omitnan')] %#ok   
MuSigFMO_NbhdInsHsp = [mean(AchievementData.minor_no_healthins_cbg(AchievementData.hisp==1),'omitnan')] %#ok   
MuSigFMO_NbhdInsWA = [mean(AchievementData.minor_no_healthins_cbg(AchievementData.black==0&AchievementData.hisp==0),'omitnan')] %#ok   
MuSigFMO_D1 = [mean(AchievementData.Dist1,'omitnan') std(AchievementData.Dist1,'omitnan') mean(isnan(AchievementData.Dist1))] %#ok
MuSigFMO_D1Blk = [mean(AchievementData.black(AchievementData.Dist1==1),'omitnan')] %#ok
MuSigFMO_D1Hsp = [mean(AchievementData.hisp(AchievementData.Dist1==1),'omitnan')] %#ok
MuSigFMO_D1WA = [1-MuSigFMO_D1Blk-MuSigFMO_D1Hsp] %#ok
MuSigFMO_D2 = [mean(AchievementData.Dist2,'omitnan') std(AchievementData.Dist2,'omitnan') mean(isnan(AchievementData.Dist2))] %#ok
MuSigFMO_D2Blk = [mean(AchievementData.black(AchievementData.Dist2==1),'omitnan')] %#ok
MuSigFMO_D2Hsp = [mean(AchievementData.hisp(AchievementData.Dist2==1),'omitnan')] %#ok
MuSigFMO_D2WA = [1-MuSigFMO_D2Blk-MuSigFMO_D2Hsp] %#ok
MuSigFMO_D3 = [mean(AchievementData.Dist3,'omitnan') std(AchievementData.Dist3,'omitnan') mean(isnan(AchievementData.Dist3))] %#ok
MuSigFMO_D3Blk = [mean(AchievementData.black(AchievementData.Dist3==1),'omitnan')] %#ok
MuSigFMO_D3Hsp = [mean(AchievementData.hisp(AchievementData.Dist3==1),'omitnan')] %#ok
MuSigFMO_D3WA = [1-MuSigFMO_D3Blk-MuSigFMO_D3Hsp] %#ok
MuSigFMO_HelpAdult = [mean(AchievementData.helpadult,'omitnan') std(AchievementData.helpadult,'omitnan') mean(isnan(AchievementData.helpadult))] %#ok
MuSigFMO_HelpPeer = [mean(AchievementData.helppeer,'omitnan') std(AchievementData.helppeer,'omitnan') mean(isnan(AchievementData.helppeer))] %#ok
MuSigFMO_NGaming = [mean(AchievementData.ngamingsys,'omitnan') std(AchievementData.ngamingsys,'omitnan') mean(isnan(AchievementData.ngamingsys))] %#ok
MuSigFMO_GamWkdy = [mean(AchievementData.gaming_weekday,'omitnan') std(AchievementData.gaming_weekday,'omitnan') mean(isnan(AchievementData.gaming_weekday))] %#ok
MuSigFMO_WkDyNet = [mean(AchievementData.weekdaynet,'omitnan') std(AchievementData.weekdaynet,'omitnan') mean(isnan(AchievementData.weekdaynet))] %#ok
MuSigFMO_Sports = [mean(AchievementData.sports,'omitnan') std(AchievementData.sports,'omitnan') mean(isnan(AchievementData.sports))] %#ok
MuSigFMO_Music = [mean(AchievementData.music,'omitnan') std(AchievementData.music,'omitnan') mean(isnan(AchievementData.music))] %#ok
MuSigFMO_Clubs = [mean(AchievementData.otheractivity,'omitnan') std(AchievementData.otheractivity,'omitnan') mean(isnan(AchievementData.otheractivity))] %#ok
MuSigFMO_FracAdSuprv = [mean((AchievementData.time_social_str)./(AchievementData.time_screen+AchievementData.time_social_unstr+AchievementData.time_social_str),'omitnan') std((AchievementData.time_social_str)./(AchievementData.time_screen+AchievementData.time_social_unstr+AchievementData.time_social_str),'omitnan') mean(isnan((AchievementData.time_social_str)./(AchievementData.time_screen+AchievementData.time_social_unstr+AchievementData.time_social_str)))] %#ok  
MuSigFMO_Intr = [mean(AchievementData.intrinsic,'omitnan') std(AchievementData.intrinsic,'omitnan') mean(isnan(AchievementData.intrinsic))] %#ok
MuSigFMO_Extr = [mean(AchievementData.extrinsic,'omitnan') std(AchievementData.extrinsic,'omitnan') mean(isnan(AchievementData.extrinsic))] %#ok
MuSigFMO_MathFav = [mean(AchievementData.math_favorite,'omitnan') std(AchievementData.math_favorite,'omitnan') mean(isnan(AchievementData.math_favorite))] %#ok
MuSigFMO_MathFavBlk = [mean(AchievementData.math_favorite(AchievementData.black==1),'omitnan')] %#ok
MuSigFMO_MathFavHsp = [mean(AchievementData.math_favorite(AchievementData.hisp==1),'omitnan')] %#ok
MuSigFMO_MathFavWA = [mean(AchievementData.math_favorite(AchievementData.black==0&AchievementData.hisp==0),'omitnan')] %#ok
MuSigFMO_MathFavURM = [mean(AchievementData.math_favorite(AchievementData.black==1|AchievementData.hisp==1),'omitnan')] %#ok
[h_MathFavDiff,p_MathFavDiff] = ttest2(AchievementData.math_favorite(AchievementData.black==1|AchievementData.hisp==1),AchievementData.math_favorite(AchievementData.black==0&AchievementData.hisp==0)) %#ok
MuSigFMO_MathLst = [mean(AchievementData.math_least,'omitnan') std(AchievementData.math_least,'omitnan') mean(isnan(AchievementData.math_least))] %#ok
MuSigFMO_STEMFav = [mean(AchievementData.STEMfav,'omitnan') std(AchievementData.STEMfav,'omitnan') mean(isnan(AchievementData.STEMfav))] %#ok
MuSigFMO_STEMFavBlk = [mean(AchievementData.STEMfav(AchievementData.black==1),'omitnan')] %#ok
MuSigFMO_STEMFavHsp = [mean(AchievementData.STEMfav(AchievementData.hisp==1),'omitnan')] %#ok
MuSigFMO_STEMFavWA = [mean(AchievementData.STEMfav(AchievementData.black==0&AchievementData.hisp==0),'omitnan')] %#ok
MuSigFMO_STEMFavURM = [mean(AchievementData.STEMfav(AchievementData.black==1|AchievementData.hisp==1),'omitnan')] %#ok
[h_STEMFavDiff,p_STEMFavDiff] = ttest2(AchievementData.STEMfav(AchievementData.black==1|AchievementData.hisp==1),AchievementData.STEMfav(AchievementData.black==0&AchievementData.hisp==0)) %#ok


disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%-----------------------         END OF TABLE 7         -------------------------------')
disp('%%%%%%--------------------------------------------------------------------------------------')
disp('%%%%%%--------------------------------------------------------------------------------------')





%%%%------------------------------------------------------------------------------------------------
%%%%------------------------------------------------------------------------------------------------
%%%%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%%%%$$$$$$$$$$$$$  SWITCH BOARD CONTROLS FOR TYPE DECOMPOSITION REGRESSIONS:  $$$$$$$$$$$$$$$$$$$$$$
%%%%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
Model            = 4.05 ;
EqualCorrelation = 0 ;
EqualVariance    = 0 ; %%%%This controls whether the variances of ThetaE and ThetaL are constrained to be the same for gender and race groups (EqualVariance==1) or not (EqualVariance==0).
ImputeMissingXs  = 1 ; %%%%This controls whether we drop rows with a missing X (ImputeMissingXs==0) or whether we impute them by E[Xmissing|{X|not missing}] (ImputeMissingXs==1).
UseGlobal        = 1 ; %%%%This controls whether we use patternsearch (UseGlobal==1) or a gradient-based algorithm in fmincon (UseGlobal==0)
%%%%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%%%%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%%%%------------------------------------------------------------------------------------------------
%%%%------------------------------------------------------------------------------------------------
NbhdVars = 2 ;
Xid              = AchievementData.id ;
inc_standardized = (AchievementData.mean_inc_cbg-mean(AchievementData.mean_inc_cbg,'omitnan'))/std(AchievementData.mean_inc_cbg,'omitnan') ;
ins_standardized = (AchievementData.minor_no_healthins_cbg) ;
ins_standardized = (ins_standardized-mean(ins_standardized,'omitnan'))/std(ins_standardized,'omitnan') ;
AchievementData.time_study_classavg = nan( size(AchievementData.time_study) ) ;


X = [inc_standardized ins_standardized ... %%%%COLUMNS 1,...,NbhdVars: Socioeconomic proxies using principal components analysis
    AchievementData.ngamingsys ... %%%%COLUMN NbhdVars+1: The number of (self-reported) gaming systems at the kid's home
    AchievementData.gaming_weekday ... %%%%COLUMN NbhdVars+2: binary for whether the kid is permitted to play videogames on weekdays
    AchievementData.weekdaynet ... %%%%COLUMN NbhdVars+3: how many daily hours of recreational internet the kid is allowed on weekdays
    AchievementData.helpadult ... %%%%COLUMN NbhdVars+4: how many adult helpers the kid has for homework
    AchievementData.helppeer ... %%%%COLUMN NbhdVars+5: how many peer helpers the kid has for homework
    AchievementData.math_favorite ... %%%%COLUMN NbhdVars+6: binary for whether math is the favorite subject
    AchievementData.math_least ... %%%%COLUMN NbhdVars+7: binary for whether math is the least favorite subject
    (AchievementData.time_social_str)./(AchievementData.time_screen+AchievementData.time_social_unstr+AchievementData.time_social_str) ... %%%%COLUMN NbhdVars+8: fraction of a student's leisure time that is under adult supervision
    AchievementData.sports  ... %%%%COLUMN NbhdVars+9: binary for whether the kid is enrolled in extra-curricular sports activities
    AchievementData.music  ... %%%%COLUMN NbhdVars+10: binary for whether the kid is enrolled in extra-curricular music activities
    AchievementData.otheractivity ... %%%%COLUMN NbhdVars+11: binary for whether the kid is enrolled in extra-curricular activities other than sports and music
    AchievementData.extrinsic  ... %%%%COLUMN NbhdVars+12: extrinsically motivated index
    AchievementData.intrinsic ... %%%%COLUMN NbhdVars+13: intrinsically motivated index
    AchievementData.fem  ... %%%%COLUMN NbhdVars+14: binary for female (omitted category is male)
    AchievementData.black  ... %%%%COLUMN NbhdVars+15: binary for black (omitted category is white/asian/other)
    AchievementData.hisp ... %%%%COLUMN NbhdVars+16: binary for hispanic (omitted category is white/asian/other)
    AchievementData.Dist2  ... %%%%COLUMN NbhdVars+17: binary for attends flossmoor school district (omitted category is wilmette)
    AchievementData.Dist3 ...  %%%%COLUMN NbhdVars+18: binary for attends chicago heights school district (omitted category is wilmette)
    AchievementData.mathpref.*AchievementData.ngamingsys ... %%%%COLUMN NbhdVars+19: Interaction between preferences for math as a subject and the consumption variables
    AchievementData.mathpref.*AchievementData.gaming_weekday ... %%%%COLUMN NbhdVars+20: Interaction between preferences for math as a subject and the consumption variables
    AchievementData.mathpref.*AchievementData.weekdaynet ... %%%%COLUMN NbhdVars+21: Interaction between preferences for math as a subject and the consumption variables
    AchievementData.grade5 ... %%%%COLUMN NbhdVars+22: Age cohort (binary for 5th grade)
    ones(length(AchievementData.id),1)... %%%%COLUMN NbhdVars+23: constant
    AchievementData.fem.*AchievementData.black... %%%%COLUMN NbhdVars+24: gender-black interaction
    AchievementData.fem.*AchievementData.hisp... %%%%COLUMN NbhdVars+25: gender-hispanic interaction
    AchievementData.fem.*AchievementData.Dist2... %%%%COLUMN NbhdVars+26: gender-flossmoor interaction
    AchievementData.fem.*AchievementData.Dist3...   %%%%COLUMN NbhdVars+27: gender-chicago heights interaction
    AchievementData.time_screen_pre...   %%%%COLUMN NbhdVars+28: daily avg screen time (recreational)
    AchievementData.mathpref.*AchievementData.time_screen_pre...   %%%%COLUMN NbhdVars+29: Interaction between preferences for math as a subject and the daily avg screen time (recreational)
    AchievementData.no_internet...   %%%%COLUMN NbhdVars+30: self-reported no stable internet connection at home
    AchievementData.time_study...   %%%%COLUMN NbhdVars+31: self-reported mean daily study time for regular schooling
    AchievementData.mobilefrac...  %%%%COLUMN NbhdVars+32: fraction of total page loads from a non-tablet mobile device (i.e., small screen)
    AchievementData.tabletfrac...  %%%%COLUMN NbhdVars+33: fraction of total page loads from a non-computer tablet device (i.e., medium size screen)
    AchievementData.understand_better...  %%%%COLUMN NbhdVars+34: This is a dummy for whether student responded "I understand [math] a lot better [relative to pre-sample period]" (omitted category is "I understand math same or less")
    AchievementData.understand_lbetter...  %%%%COLUMN NbhdVars+35:This is a dummy for whether student responded "I understand [math] a little better [relative to pre-sample period]" (omitted category is "I understand math same or less")
    AchievementData.enjoy_better...  %%%%COLUMN NbhdVars+36: This is a dummy for whether the student responded "I like math a lot more now [relative to pre-sample period]"  (omitted category is "I enjoy math same or less") 
    AchievementData.enjoy_lbetter...  %%%%COLUMN NbhdVars+37: This is a dummy for whether the student responded "I like math a little better now [relative to pre-sample period]"  (omitted category is "I enjoy math same or less")
    AchievementData.miss_hmwk_some...  %%%%COLUMN NbhdVars+38: This is a dummy for whether the student reported missing homework sometimes but less than half the time   (omitted category is "I never miss homework assignments")
    AchievementData.miss_hmwk_usually...  %%%%COLUMN NbhdVars+39: This is a dummy for whether the student reported missing homework sometimes half or more of the time   (omitted category is "I never miss homework assignments")
    AchievementData.involved_parent...  %%%%COLUMN NbhdVars+40: <<FROM PARENT SURVEY---333OBS>> dummy for whether parent self-reported positive time spent daily with child on schoolwork exceeds 2hrs on average
    AchievementData.bigfam...  %%%%COLUMN NbhdVars+41: <<FROM PARENT SURVEY---333OBS>> dummy for # of children in the household>=3
    AchievementData.middle_child...  %%%%COLUMN NbhdVars+42: <<FROM PARENT SURVEY---307OBS>> This is a dummy for middle child status at the time of the sample period ("oldest" means oldest OR only child)
    AchievementData.youngest_child] ;  %%%%COLUMN NbhdVars+43: <<FROM PARENT SURVEY---307OBS>> This is a dummy for youngest child status at the time of the sample period ("oldest" means oldest OR only child)
clear wtsL wtsE


Y_E = AchievementData.ThetaE ;
Y_L = AchievementData.ThetaL ;
Q   = AchievementData.Q ;
con = AchievementData.contract ;
fem = AchievementData.fem ;
blk = AchievementData.black ;
hsp = AchievementData.hisp ;

if ImputeMissingXs==1  %%%%Here we will impute missing X values by regressing the missing X on the other Xs, using the subsample where both are available.
    disp('Running Imputation Algorithm for missing Xs...')
    Xtemp = X ;  %%%%We make a second copy because while we are doing imputation, we want to be able to tell which rows were complete in the original dataset.
    allvars = (1:length(Xtemp(1,:))) ;
    for ii=1:length(Xtemp(:,1))
        missingvars = find(isnan(Xtemp(ii,:))) ;
        if ~isempty(missingvars)
            availablevars = setdiff(allvars,missingvars) ;
            for jj=1:length(missingvars)
                Ytemp2 = X(:,missingvars(jj)) ;
                Xtemp2 = X(:,availablevars) ;
                dropindx = find(isnan(sum([Ytemp2,Xtemp2],2))) ;
                Ytemp2(dropindx) = [] ;
                Xtemp2(dropindx,:) = [] ;
                warning('error', 'MATLAB:singularMatrix');
                Bimpute = (Xtemp2'*Xtemp2)\Xtemp2'*Ytemp2 ;  %regress(Ytemp2,Xtemp2) ;
                Xtemp(ii,missingvars(jj)) = Xtemp(ii,availablevars)*Bimpute ;
            end
        end
    end
    X = Xtemp ;
    useindx  = find(~isnan(sum(X,2))) ;
    missindx = find(isnan(sum(X,2))) ;
    clear Xtemp allvars missingvars available vars ii jj Ytemp2 Xtemp2 dropindx Bimpute
else
    useindx  = find(~isnan(sum(X,2))) ;
    missindx = find(isnan(sum(X,2))) ;
end
frac_no_missing = 1-length(missindx)/(length(useindx)+length(missindx)) %#ok

if Model==1   %%%%SPECIFICATION 1 IN TABLES 2 AND 3
              %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Grade-5, 7 CONST
    colindx_E = [1:NbhdVars,NbhdVars+14:NbhdVars+16,NbhdVars+22:NbhdVars+23] ;
    X_E       = X(:,colindx_E) ;
    X_L       = X(:,colindx_E) ;
    clear colindx_E
elseif Model==2  %%%%SPECIFICATION 2 IN TABLES 2 AND 3
                 %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
    colindx_E = [1:NbhdVars,NbhdVars+14:NbhdVars+14+4,NbhdVars+22:NbhdVars+23] ;
    X_E       = X(:,colindx_E) ;
    X_L       = X(:,colindx_E) ;
    clear colindx_E
elseif Model==3 %%%%%SPECIFICATION 3 IN TABLES 2 AND 3
                %%%%theta_E Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                %%%%                   5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5, 11 CONST
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                %%%%                    5 MathFav, 6 MathLeast, 7 EXTRINSIC, 8 INTRINSIC, 9 fem, 
                %%%%                    10 Black, 11 Hispanic, 12 Dist2, 13 Dist3, 14 Grade-5, 15 CONST
    colindx_E = [1:NbhdVars,NbhdVars+4:NbhdVars+4+1,NbhdVars+14:NbhdVars+14+4,NbhdVars+22:NbhdVars+23] ;
    X_E       = X(:,colindx_E) ;
    colindx_L = [1:NbhdVars,NbhdVars+4:NbhdVars+7,NbhdVars+12:NbhdVars+18,NbhdVars+22:NbhdVars+23] ;
    X_L       = X(:,colindx_L) ;
    clear colindx_E
elseif Model==4.05 %%%%%SPECIFICATION 4 IN TABLES 2 AND 3
                   %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                   %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                   %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                   %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
                   %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
                   %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                   %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
                   %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                   %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
    colindx_E = [1:NbhdVars,NbhdVars+4:NbhdVars+4+1,NbhdVars+14:NbhdVars+14+4,NbhdVars+22:NbhdVars+23,NbhdVars+30,NbhdVars+32:NbhdVars+33] ;
    X_E       = X(:,colindx_E) ;
    X_L       = X(:,[1:NbhdVars+18,NbhdVars+22:NbhdVars+23,NbhdVars+31,NbhdVars+28,NbhdVars+30,NbhdVars+32:NbhdVars+33]) ;  %%%%Here we use all variables but the interactions between Math preferences and the gaming/internet consumption variables (in positions NbhdVars+19:NbhdVars+21).
    clear colindx_E
elseif Model==4.085%%%%%SPECIFICATION 5 IN TABLES 2 AND 3
                   %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                   %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                   %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                   %%%%                    15 involved_par, 16 bigfam, 17 middlechild, 18 youngestchild
                   %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
                   %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
                   %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                   %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
                   %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                   %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
                   %%%%                    28 involved_par, 29 bigfam, 30 middlechild, 31 youngestchild
    colindx_E = [1:NbhdVars,NbhdVars+4:NbhdVars+4+1,NbhdVars+14:NbhdVars+14+4,NbhdVars+22:NbhdVars+23,NbhdVars+30,NbhdVars+32:NbhdVars+33,NbhdVars+40,NbhdVars+41,NbhdVars+42:NbhdVars+43] ; 
    X_E       = X(:,colindx_E) ; 
    X_L       = X(:,[1:NbhdVars+18,NbhdVars+22:NbhdVars+23,NbhdVars+31,NbhdVars+28,NbhdVars+30,NbhdVars+32:NbhdVars+33,NbhdVars+40,NbhdVars+41,NbhdVars+42:NbhdVars+43]) ;  
    clear colindx_E
end




FIRST_PASS = 0 ;
if Model==1
    % load('betareg0_Model1_Impute1.mat')
    load('.\STAGE3\betareg0_Model1_Impute1Restart.mat')
    %%%%This is the basic version of the regression model: 
    %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Grade-5, 7 CONST
elseif Model==2
    if FIRST_PASS==1
        load('.\STAGE3\betareg0_Model1_Impute1Restart.mat')
        betareg0 = [betareg0(1:5,:);
            zeros(2,length(betareg0(1,:)));
            betareg0(6:12,:);
            zeros(2,length(betareg0(1,:)));
            betareg0(13:end,:)] ;
    else
        load('.\STAGE3\betareg0_Model2_Impute1Restart.mat')
    end
    %%%%This is the basic version of the regression model: 
    %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
elseif Model==3
    if FIRST_PASS==1
        load('.\STAGE3\betareg0_Model2_Impute1Restart.mat')
        betareg0 = [betareg0(1:2,:);
            zeros(2,length(betareg0(1,:)));
            betareg0(3:11,:);
            zeros(6,length(betareg0(1,:)));
            betareg0(12:end,:)] ;
    else
        load('.\STAGE3\betareg0_Model3_Impute1Restart.mat')
    end
    %%%%This is the basic version of the regression model: 
    %%%%theta_E Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                %%%%                   5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5, 11 CONST
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                %%%%                    5 MathFav, 6 MathLeast, 7 EXTRINSIC, 8 INTRINSIC, 9 fem, 
                %%%%                    10 Black, 11 Hispanic, 12 Dist2, 13 Dist3, 14 Grade-5, 15 CONST
elseif Model==4.05   
    % load('betareg0_Model4point05_Impute1.mat') 
    % betareg0 = [betareg0(1:12,1);0;0;betareg0(13:35,1);0;0;betareg0(36:end,1)] ;
    load('.\STAGE3\betareg0_Model4point05_Impute1Restart.mat')
    %%%%This is the basic version of the regression model: 
                   %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                   %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                   %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                   %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
                   %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
                   %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                   %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
                   %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                   %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
elseif Model==4.085   
    if FIRST_PASS==1
        load('.\STAGE3\betareg0_Model4point05_Impute1Restart.mat')
        betareg0 = [betareg0(1:14,:);
            zeros(5,length(betareg0(1,:)));
            betareg0(15:41,:);
            zeros(5,length(betareg0(1,:)));
            betareg0(42:end,:)] ;
    else
        load('.\STAGE3\betareg0_Model4point085_Impute1Restart.mat')
    end
    %%%%This is the basic version of the regression model: 
    %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
    %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
    %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
    %%%%                    15 Involved_par, 16 bigfam, 17 middlechild, 18 youngestchild
    %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
    %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
    %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
    %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
    %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
    %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
    %%%%                    28 Involved_par, 29 bigfam, 30 middlechild, 31 youngestchild
end




%%%%% FOR THE SOLUTION WE USE A COMBINATION OF QUASI-NEWTON AND DERIVATIVE-FREE METHODS.  WE ALSO
%%%%% USE A MULTIPLE RE-START APPROACH TO ENSURE A RELIABLE SOLUTION.
RunTobitEstimator = 0 ;
if RunTobitEstimator==1
    diary('StressTestTobitSolution.txt')  


    tic
    for ll=1:1 %length(betareg0(1,:))
        [betaregtemp,MLEouttemp] = ...
            TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
            Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
            LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
            betareg0(:,ll),wts(useindx),0.1,EqualVariance,EqualCorrelation,1,10^-8) ;  %%%%Interior-point algorithm
        logLtempvals = [MLEouttemp.logLval] %#ok
        [betaregtemp,MLEouttemp] = ...
            TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
            Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
            LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
            betaregtemp,wts(useindx),0.2,EqualVariance,EqualCorrelation,1,10^-8) ;  %%%%SQP-Legacy algorithm
        logLtempvals = [MLEouttemp.logLval logLtempvals] %#ok
        [betaregtemp,MLEouttemp] = ...
            TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
            Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
            LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
            betaregtemp,wts(useindx),0.3,EqualVariance,EqualCorrelation,1,10^-8) ;  %%%%SQP algorithm
        logLtempvals = [MLEouttemp.logLval logLtempvals] %#ok
        [betaregtemp,MLEouttemp] = ...
            TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
            Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
            LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
            betaregtemp,wts(useindx),0.4,EqualVariance,EqualCorrelation,1,10^-8) ;  %%%%Active-Set algorithm
        logLtempvals = [MLEouttemp.logLval logLtempvals] %#ok
        betareg0(:,ll) = betaregtemp ;
        if ll==1
            MLEout  = MLEouttemp ; %#ok
            betareg = betaregtemp ; %#ok
        end
    end

    [betareg4,MLEout4] = ...
        TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
        Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
        LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
        betareg0(:,1),wts(useindx),4,EqualVariance,EqualCorrelation,1,10^-6) ;  %%%%Patternsearch with NUPS-MADS
    [betareg1,MLEout1] = ...
        TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
        Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
        LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
        betareg4,wts(useindx),1,EqualVariance,EqualCorrelation,1,10^-6) ;  %%%%Patternsearch with CLASSIC
    [betareg2,MLEout2] = ...
        TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
        Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
        LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
        betareg1,wts(useindx),2,EqualVariance,EqualCorrelation,1,10^-6) ;  %%%%Patternsearch with NUPS
    [betareg3,MLEout3] = ...
        TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
        Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
        LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
        betareg2,wts(useindx),3,EqualVariance,EqualCorrelation,1,10^-6) ;  %%%%Patternsearch with NUPS-GPS
    [betareg,MLEout] = ...
        TypeRegressionEstimatorSmoothed(Y_E(useindx),X_E(useindx,:),Y_L(useindx),X_L(useindx,:),...
        Q(useindx),con(useindx),fem(useindx),blk(useindx),hsp(useindx),...
        LPointEst.domE,LPointEst.threshold1,LPointEst.threshold2,LPointEst.threshold3,...
        betareg3,wts(useindx),0.1,EqualVariance,EqualCorrelation,1,10^-12) ;
    betareg0(:,1) = betareg ;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%THIS NEXT LINE OF CODE EXECUTES THE MULTI-METHOD MULTIPLE RESTART SOLUTION PROCESS
    %%%%%%%%%%%DESCRIBED IN APPENDIX A.5
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%Here we test our solution against multiple restarts:
    run StressTestTobitSolution
    userInput = input('YOU HAVE NOW COMPLETED A ROUND OF MULTIPLE RESTARTS FOR THE TOBIT ESTIMATOR.  IF THE CURRENT ROUND RESULTED IN A SIGNIFICANT IMPROVEMENT ON THE PREVIOUS START POINT (FIRST COLUMN OF betareg0), THEN ENTER 1.  IF YOU DO SO, COMPUTATION WILL TERMINATE, AND YOU CAN HIT "Run Section" IN THE MATLAB EDITOR TO START A NEW ROUND.  OTHERWISE, ENTER 0, AND COMPUTATION WILL CONTINUE TO THE NEXT STAGE.')
    if userInput==1
        error('TobitComputation:OneMoreRound','Computation halted temporarily.  Hit "Run Section" in the MATLAB editor to re-execute the block of Tobit estimation code beginning on line 904 of AAA_File00_Driver.m.')
    elseif userInput==0
        disp('Now saving Tobit Estimation Results...')
    else
        error('TobitComputation:OneMoreRound','You must enter either 1 or 0...')
    end
    diary off


    if Model==1
        save('.\STAGE3\betareg0_Model1_Impute1Restart.mat','betareg0','betareg','MLEout')
    elseif Model==2
        save('.\STAGE3\betareg0_Model2_Impute1Restart.mat','betareg0','betareg','MLEout')
    elseif Model==3
        save('.\STAGE3\betareg0_Model3_Impute1Restart.mat','betareg0','betareg','MLEout')
    elseif Model==4.05
        save('.\STAGE3\betareg0_Model4point05_Impute1Restart.mat','betareg0','betareg','MLEout')
    elseif Model==4.085
        save('.\STAGE3\betareg0_Model4point085_Impute1Restart.mat','betareg0','betareg','MLEout')
    end
else
    if Model==1
        load('.\STAGE3\betareg0_Model1_Impute1Restart.mat')
    elseif Model==2
        load('.\STAGE3\betareg0_Model2_Impute1Restart.mat')
    elseif Model==3
        load('.\STAGE3\betareg0_Model3_Impute1Restart.mat')
    elseif Model==4.05
        load('.\STAGE3\betareg0_Model4point05_Impute1Restart.mat')
    elseif Model==4.085
        load('.\STAGE3\betareg0_Model4point085_Impute1Restart.mat')
    end
end

Vhat   = inv(MLEout.hess) ;  
StdErr = sqrt(diag(Vhat)) ;

disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%           TABLES 2 & 3 WALD TESTS:        %%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')

%%%%THE FOLLOWING SECTION OF CODE RUNS A TEST FOR EXCLUSION RESTRICTIONS:
RunWaldTest = 1 ;
if RunWaldTest==1
    disp('-----------------------------------------------------------------------------------------------')
    disp('-----------------------------------------------------------------------------------------------')
    disp('>>')
    disp('RUNNING STAGE-3 WALD TESTS FOR JOINT EXCLUSIONS FROM THE TOBIT MODEL...')
    disp('>>')
    disp('-----------------------------------------------------------------------------------------------')
    disp('-----------------------------------------------------------------------------------------------')
    disp('>>')
    logL_u = -1*MLEout.logLval ;
    for tt=1:5 %%%run 3 tests: one each for exclusion of school FEs from the log(thetaE) and
        %%%%log(ThetaL) equations individually, and a third for joint exclusion of school
        %%%%FEs from both equations at once.  The 4th/5th rounds are a battery of Wald tests
        %%%%for variables appearing in the yes/no sections of tables 2 and 3
        if tt==1 && Model>=2
            disp('>>')
            disp('>>')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('$$$$$$$$  TEST 1: JOINT EXCLUSION OF SCHOOL FEs FROM log(thetaE) EQUATION...  $$$$$$$')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('>>')
            if Model==2
                %%%%This is the basic version of the regression model:
                %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
                R = [zeros(1,5), 1, 0, 0, 0, zeros(1,9), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,5), 0, 1, 0, 0, zeros(1,9), zeros(1,12)]; %%%%zero out Dist3 from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==3
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                   5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5, 11 CONST
                %%%%theta_L Variables: 15
                R = [zeros(1,7), 1, 0, 0, 0, zeros(1,15), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, 0, 0, zeros(1,15), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==4.05
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                %%%%theta_L Variables: 27
                R = [zeros(1,7), 1, 0, zeros(1,5), zeros(1,27), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, zeros(1,5), zeros(1,27), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==4.085
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                %%%%                    15 Involved_par, 16 bigfam, 17 middlechild, 18 youngestchild
                %%%%theta_L Variables: 31
                R = [zeros(1,7), 1, 0, zeros(1,9), zeros(1,31), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, zeros(1,9), zeros(1,31), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaE)
                r = zeros(2,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
        elseif tt==2 && Model>=2
            disp('>>')
            disp('>>')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('$$$$$$$$  TEST 2: JOINT EXCLUSION OF SCHOOL FEs FROM log(thetaL) EQUATION...  $$$$$$$')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('>>')
            if Model==2
                %%%%This is the basic version of the regression model:
                %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
                R = [zeros(1,9), zeros(1,5), 1, 0, 0, 0, zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,9), zeros(1,5), 0, 1, 0, 0, zeros(1,12)]; %%%%zero out Dist3 from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==3
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables: 11
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 MathFav, 6 MathLeast, 7 EXTRINSIC, 8 INTRINSIC, 9 fem,
                %%%%                    10 Black, 11 Hispanic, 12 Dist2, 13 Dist3, 14 Grade-5, 15 CONST
                R = [zeros(1,11), zeros(1,11), 1, 0, 0, 0, zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,11), zeros(1,11), 0, 1, 0, 0, zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.05
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  14
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay,
                %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav,
                %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem,
                %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
                R = [zeros(1,14), zeros(1,18), 1, 0, zeros(1,7), zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,14), zeros(1,18), 0, 1, zeros(1,7), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.085
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  18
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay,
                %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav,
                %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem,
                %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
                %%%%                    28 Involved_par, 29 bigfam, 30 middlechild, 31 youngestchild
                R = [zeros(1,18), zeros(1,18), 1, 0, zeros(1,11), zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,18), zeros(1,18), 0, 1, zeros(1,11), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(2,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
        elseif tt==3 && Model>=2
            disp('>>')
            disp('>>')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('$$$$$$$$$$  TEST 3: JOINT EXCLUSION OF SCHOOL FEs FROM BOTH EQUATIONS...  $$$$$$$$$$$')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('>>')
            if Model==2
                %%%%This is the basic version of the regression model:
                %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
                R = [zeros(1,5), 1, 0, 0, 0, zeros(1,9), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,5), 0, 1, 0, 0, zeros(1,9), zeros(1,12);... %%%%zero out Dist3 from log(thetaE)
                    zeros(1,9), zeros(1,5), 1, 0, 0, 0, zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,9), zeros(1,5), 0, 1, 0, 0, zeros(1,12)]; %%%%zero out Dist3 from log(thetaL)
                r = zeros(4,1) ;
            elseif Model==3
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                   5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5, 11 CONST
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 MathFav, 6 MathLeast, 7 EXTRINSIC, 8 INTRINSIC, 9 fem,
                %%%%                    10 Black, 11 Hispanic, 12 Dist2, 13 Dist3, 14 Grade-5, 15 CONST
                R = [zeros(1,7), 1, 0, 0, 0, zeros(1,15), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, 0, 0, zeros(1,15), zeros(1,12);... %%%%zero out Dist3 from log(thetaE)
                    zeros(1,11), zeros(1,11), 1, 0, 0, 0, zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,11), zeros(1,11), 0, 1, 0, 0, zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(4,1) ;
            elseif Model==4.05
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay,
                %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav,
                %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem,
                %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
                R = [zeros(1,7), 1, 0, zeros(1,5), zeros(1,27), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, zeros(1,5), zeros(1,27), zeros(1,12);... %%%%zero out Dist3 from log(thetaE)
                    zeros(1,14), zeros(1,18), 1, 0, zeros(1,7), zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,14), zeros(1,18), 0, 1, zeros(1,7), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(4,1) ;
            elseif Model==4.085
                %%%%This is the basic version of the regression model:
                %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER,
                %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
                %%%%                    15 Involved_par, 16 bigfam, 17 middlechild, 18 youngestchild
                %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay,
                %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav,
                %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem,
                %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
                %%%%                    28 Involved_par, 29 bigfam, 30 middlechild, 31 youngestchild
                R = [zeros(1,7), 1, 0, zeros(1,9), zeros(1,31), zeros(1,12);... %%%%zero out Dist2 from log(thetaE)
                    zeros(1,7), 0, 1, zeros(1,9), zeros(1,31), zeros(1,12);... %%%%zero out Dist3 from log(thetaE)
                    zeros(1,18), zeros(1,18), 1, 0, zeros(1,11), zeros(1,12);... %%%%zero out Dist2 from log(thetaL)
                    zeros(1,18), zeros(1,18), 0, 1, zeros(1,11), zeros(1,12)] ; %%%%zero out Dist3 from log(thetaL)
                r = zeros(4,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
        elseif tt==4
            disp('>>')
            disp('>>')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('$$$$$$$$$$$$  TEST 4: JOINT EXCLUSIONS OF "YES/NO" VARIABLE GROUPS FROM log(thetaE)...  $$$$$$$$$$$$$$')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('>>')

            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('EXCLUSION OF NBHD-SES CONTROLS:')
            disp('------------------------------------------------------------------------')
            disp('>>')
            if Model==1
                R = [1, 0, zeros(1,5), zeros(1,7), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaE)
                    0, 1, zeros(1,5), zeros(1,7), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==2
                R = [1, 0, zeros(1,7), zeros(1,9), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaE)
                    0, 1, zeros(1,7), zeros(1,9), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==3
                R = [1, 0, zeros(1,9), zeros(1,15), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaE)
                    0, 1, zeros(1,9), zeros(1,15), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==4.05
                R = [1, 0, zeros(1,12), zeros(1,27), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaE)
                    0, 1, zeros(1,12), zeros(1,27), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaE)
                r = zeros(2,1) ;
            elseif Model==4.085
                R = [1, 0, zeros(1,16), zeros(1,31), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaE)
                    0, 1, zeros(1,16), zeros(1,31), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaE)
                r = zeros(2,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('------------------------------------------------------------------------')
            disp('>>')

            if Model>=3
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF HOME ACADEMIC SUPPORT CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==3
                    R = [zeros(1,2), 1, 0, zeros(1,7), zeros(1,15), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaE)
                        zeros(1,2), 0, 1, zeros(1,7), zeros(1,15), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaE)
                    r = zeros(2,1) ;
                elseif Model==4.05
                    R = [zeros(1,2), 1, 0, zeros(1,10), zeros(1,27), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaE)
                        zeros(1,2), 0, 1, zeros(1,10), zeros(1,27), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaE)
                    r = zeros(2,1) ;
                elseif Model==4.085
                    R = [zeros(1,2), 1, 0, zeros(1,14), zeros(1,31), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaE)
                        zeros(1,2), 0, 1, zeros(1,14), zeros(1,31), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaE)
                    r = zeros(2,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end

            if Model>=4.05
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF HOME CONNECTIVITY CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==4.05
                    R = [zeros(1,11), 1, 0, 0, zeros(1,27), zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaE)
                        zeros(1,11), 0, 1, 0, zeros(1,27), zeros(1,12);...   %%%%zero out MOBILEFRAC from log(thetaE)
                        zeros(1,11), 0, 0, 1, zeros(1,27), zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaE)
                    r = zeros(3,1) ;
                elseif Model==4.085
                    R = [zeros(1,11), 1, 0, 0, zeros(1,4), zeros(1,31), zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaE)
                        zeros(1,11), 0, 1, 0, zeros(1,4), zeros(1,31), zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaE)
                        zeros(1,11), 0, 0, 1, zeros(1,4), zeros(1,31), zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaE)
                    r = zeros(3,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end

            if Model==4.085
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF PARENT SURVEY CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                R = [zeros(1,14), 1, 0, 0, 0, zeros(1,31), zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaE)
                    zeros(1,14), 0, 1, 0, 0, zeros(1,31), zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaE)
                    zeros(1,14), 0, 0, 1, 0, zeros(1,31), zeros(1,12);... %%%%zero out TABLETFRAC from log(thetaE)
                    zeros(1,14), 0, 0, 0, 1, zeros(1,31), zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaE)
                r = zeros(4,1) ;
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end
        elseif tt==5
            disp('>>')
            disp('>>')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('$$$$$$$$$$$$  TEST 5: JOINT EXCLUSIONS OF "YES/NO" VARIABLE GROUPS FROM log(thetaL)...  $$$$$$$$$$$$$$')
            disp('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
            disp('>>')

            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('EXCLUSION OF RACE CONTROLS:')
            disp('------------------------------------------------------------------------')
            disp('>>')
            if Model==1
                R = [zeros(1,7), 0, 0, 0, 1, 0, 0, 0, zeros(1,12);... %%%%zero out Black from log(thetaL)
                    zeros(1,7), 0, 0, 0, 0, 1, 0, 0, zeros(1,12)];   %%%%zero out Hispanic from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==2
                R = [zeros(1,9), zeros(1,3), 1, 0, zeros(1,4), zeros(1,12);... %%%%zero out Black from log(thetaL)
                    zeros(1,9), zeros(1,3), 0, 1, zeros(1,4), zeros(1,12)];   %%%%zero out Hispanic from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==3
                R = [zeros(1,11), zeros(1,9), 1, 0, zeros(1,4), zeros(1,12);... %%%%zero out Black from log(thetaL)
                    zeros(1,11), zeros(1,9), 0, 1, zeros(1,4), zeros(1,12)];   %%%%zero out Hispanic from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.05
                R = [zeros(1,14), zeros(1,16), 1, 0, zeros(1,9), zeros(1,12);... %%%%zero out Black from log(thetaL)
                    zeros(1,14), zeros(1,16), 0, 1, zeros(1,9), zeros(1,12)];   %%%%zero out Hispanic from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.085
                R = [zeros(1,18), zeros(1,16), 1, 0, zeros(1,13), zeros(1,12);... %%%%zero out Black from log(thetaL)
                    zeros(1,18), zeros(1,16), 0, 1, zeros(1,13), zeros(1,12)];   %%%%zero out Hispanic from log(thetaL)
                r = zeros(2,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('------------------------------------------------------------------------')
            disp('>>')

            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('EXCLUSION OF NBHD-SES CONTROLS:')
            disp('------------------------------------------------------------------------')
            disp('>>')
            if Model==1
                R = [zeros(1,7), 1, 0, zeros(1,5), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaL)
                    zeros(1,7), 0, 1, zeros(1,5), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==2
                R = [zeros(1,9), 1, 0, zeros(1,7), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaL)
                    zeros(1,9), 0, 1, zeros(1,7), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==3
                R = [zeros(1,11), 1, 0, zeros(1,13), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaL)
                    zeros(1,11), 0, 1, zeros(1,13), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.05
                R = [zeros(1,14), 1, 0, zeros(1,25), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaL)
                    zeros(1,14), 0, 1, zeros(1,25), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaL)
                r = zeros(2,1) ;
            elseif Model==4.085
                R = [zeros(1,18), 1, 0, zeros(1,29), zeros(1,12);... %%%%zero out NBHD-Inc from log(thetaL)
                    zeros(1,18), 0, 1, zeros(1,29), zeros(1,12)];   %%%%zero out NBHD UninsMinors from log(thetaL)
                r = zeros(2,1) ;
            end
            N        = length(AchievementData.Q) ;
            disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
            WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
            DF       = length(r) ;
            disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
            LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
            disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
            disp('>>')
            disp('------------------------------------------------------------------------')
            disp('------------------------------------------------------------------------')
            disp('>>')

            if Model>=3
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF HOME ACADEMIC SUPPORT CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==3
                    R = [zeros(1,11), zeros(1,2), 1, 0, zeros(1,11), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaL)
                        zeros(1,11), zeros(1,2), 0, 1, zeros(1,11), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.05
                    R = [zeros(1,14), zeros(1,5), 1, 0, zeros(1,20), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaL)
                        zeros(1,14), zeros(1,5), 0, 1, zeros(1,20), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,5), 1, 0, zeros(1,24), zeros(1,12);... %%%%zero out HELPER_ADULT from log(thetaL)
                        zeros(1,18), zeros(1,5), 0, 1, zeros(1,24), zeros(1,12)];   %%%%zero out HELPER_PEER from log(thetaL)
                    r = zeros(2,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end

            if Model>=4.05
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF EXTRACURRICULARS CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==4.05
                    R = [zeros(1,14), zeros(1,9), 1, 0, 0, 0, zeros(1,14), zeros(1,12);... %%%%zero out FracSocialSuperv from log(thetaL)
                        zeros(1,14), zeros(1,9), 0, 1, 0, 0, zeros(1,14), zeros(1,12);... %%%%zero out Sports from log(thetaL)
                        zeros(1,14), zeros(1,9), 0, 0, 1, 0, zeros(1,14), zeros(1,12);... %%%%zero out music from log(thetaL)
                        zeros(1,14), zeros(1,9), 0, 0, 0, 1, zeros(1,14), zeros(1,12)];   %%%%zero out OtherActiv from log(thetaL)
                    r = zeros(4,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,9), 1, 0, 0, 0, zeros(1,18), zeros(1,12);... %%%%zero out FracSocialSuperv from log(thetaL)
                        zeros(1,18), zeros(1,9), 0, 1, 0, 0, zeros(1,18), zeros(1,12);... %%%%zero out Sports from log(thetaL)
                        zeros(1,18), zeros(1,9), 0, 0, 1, 0, zeros(1,18), zeros(1,12);... %%%%zero out music from log(thetaL)
                        zeros(1,18), zeros(1,9), 0, 0, 0, 1, zeros(1,18), zeros(1,12)];   %%%%zero out OtherActiv from log(thetaL)
                    r = zeros(4,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')

                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF GAMING/SURFING CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==4.05
                    R = [zeros(1,14), zeros(1,2), 1, 0, 0, zeros(1,18), 0, zeros(1,3), zeros(1,12);... %%%%zero out #GamingSys from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 1, 0, zeros(1,18), 0, zeros(1,3), zeros(1,12);... %%%%zero out GamingWkDay from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 1, zeros(1,18), 0, zeros(1,3), zeros(1,12);... %%%%zero out WkDayNet from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,18), 1, zeros(1,3), zeros(1,12)];   %%%%zero out TimeScreen from log(thetaL)
                    r = zeros(4,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,2), 1, 0, 0, zeros(1,18), 0, zeros(1,7), zeros(1,12);... %%%%zero out #GamingSys from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 1, 0, zeros(1,18), 0, zeros(1,7), zeros(1,12);... %%%%zero out GamingWkDay from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 1, zeros(1,18), 0, zeros(1,7), zeros(1,12);... %%%%zero out WkDayNet from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,18), 1, zeros(1,7), zeros(1,12)];   %%%%zero out TimeScreen from log(thetaL)
                    r = zeros(4,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')

                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF HOME CONNECTIVITY CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==4.05
                    R = [zeros(1,14), zeros(1,24), 1, 0, 0, zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaL)
                        zeros(1,14), zeros(1,24), 0, 1, 0, zeros(1,12);...   %%%%zero out MOBILEFRAC from log(thetaL)
                        zeros(1,14), zeros(1,24), 0, 0, 1, zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaL)
                    r = zeros(3,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,24), 1, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaL)
                        zeros(1,18), zeros(1,24), 0, 1, 0, zeros(1,4), zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaL)
                        zeros(1,18), zeros(1,24), 0, 0, 1, zeros(1,4), zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaL)
                    r = zeros(3,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')

                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('JOINT EXCLUSION OF OUTSIDE TIME-USE CONTROLS: GAMING/SURFING (3), EXTRACURRICULARS (4), REG STUDY TIME, REG SCREEN TIME, HOME CONNECTIVITY (3)')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==4.05
                    R = [zeros(1,14), zeros(1,2), 1, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out #GamingSys from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 1, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out GamingWkDay from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 1, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out WkDayNet from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 1, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out FracSocialSuperv from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 1, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out Sports from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 1, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out music from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 1, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,12);... %%%%zero out OtherActiv from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 1, 0, 0, 0, 0, zeros(1,12);... %%%%zero out TimeStudy from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 1, 0, 0, 0, zeros(1,12);... %%%%zero out TimeScreen from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 1, 0, 0, zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 1, 0, zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaL)
                        zeros(1,14), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 1, zeros(1,12)]; %%%%zero out TABLETFRAC from log(thetaL)
                    r = zeros(12,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,2), 1, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out #GamingSys from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 1, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out GamingWkDay from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 1, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out WkDayNet from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 1, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out FracSocialSuperv from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 1, 0, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out Sports from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 1, 0, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out music from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 1, zeros(1,9), 0, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out OtherActiv from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 1, 0, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out TimeStudy from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 1, 0, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out TimeScreen from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 1, 0, 0, zeros(1,4), zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 1, 0, zeros(1,4), zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaL)
                        zeros(1,18), zeros(1,2), 0, 0, 0, zeros(1,4), 0, 0, 0, 0, zeros(1,9), 0, 0, 0, 0, 1, zeros(1,4), zeros(1,12)]; %%%%zero out TABLETFRAC from log(thetaL)
                    r = zeros(12,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end

            if Model==4.085
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF PARENT SURVEY CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                R = [zeros(1,18), zeros(1,27), 1, 0, 0, 0, zeros(1,12);... %%%%zero out NO_INTERNET from log(thetaL)
                    zeros(1,18), zeros(1,27), 0, 1, 0, 0, zeros(1,12);... %%%%zero out MOBILEFRAC from log(thetaL)
                    zeros(1,18), zeros(1,27), 0, 0, 1, 0, zeros(1,12);... %%%%zero out TABLETFRAC from log(thetaL)
                    zeros(1,18), zeros(1,27), 0, 0, 0, 1, zeros(1,12)];   %%%%zero out TABLETFRAC from log(thetaL)
                r = zeros(4,1) ;
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end

            if Model>=3
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF SUBJECT PREFERENCE CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==3
                    R = [zeros(1,11), zeros(1,4), 1, 0, zeros(1,9), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,11), zeros(1,4), 0, 1, zeros(1,9), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.05
                    R = [zeros(1,14), zeros(1,7), 1, 0, zeros(1,18), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,14), zeros(1,7), 0, 1, zeros(1,18), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,7), 1, 0, zeros(1,22), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,18), zeros(1,7), 0, 1, zeros(1,22), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')


                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('EXCLUSION OF MOTIVATION TYPE CONTROLS:')
                disp('------------------------------------------------------------------------')
                disp('>>')
                if Model==3
                    R = [zeros(1,11), zeros(1,6), 1, 0, zeros(1,7), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,11), zeros(1,6), 0, 1, zeros(1,7), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.05
                    R = [zeros(1,14), zeros(1,13), 1, 0, zeros(1,12), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,14), zeros(1,13), 0, 1, zeros(1,12), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                elseif Model==4.085
                    R = [zeros(1,18), zeros(1,13), 1, 0, zeros(1,16), zeros(1,12);... %%%%zero out MathFav from log(thetaL)
                        zeros(1,18), zeros(1,13), 0, 1, zeros(1,16), zeros(1,12)];   %%%%zero out MathLeast from log(thetaL)
                    r = zeros(2,1) ;
                end
                N        = length(AchievementData.Q) ;
                disp(strcat('MODEL=',num2str(Model),'... and sample size is N=',num2str(N)))
                WaldStat = (R*betareg - r)'*inv(R*Vhat*R')*(R*betareg - r) ; %#ok
                DF       = length(r) ;
                disp(strcat('WaldStat=',num2str(WaldStat),'... and degrees of freedom are DF=',num2str(DF)))
                LR_Pval  = 1 - chi2cdf(WaldStat,DF) ;
                disp(strcat('TEST RESULTS IN A P-VALUE OF>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>',num2str(LR_Pval),'.'))
                disp('>>')
                disp('------------------------------------------------------------------------')
                disp('------------------------------------------------------------------------')
                disp('>>')
            end
        end

    end
    disp('-----------------------------------------------------------------------------------------------')
    disp('-----------------------------------------------------------------------------------------------')
    disp('-----------------------------------------------------------------------------------------------')
    disp('-----------------------------------------------------------------------------------------------')
end
    




 

           
format shortg                      

MLEout  %#ok
endbetas   = length(X_E(1,:)) + length(X_L(1,:)) ;  %%%%This is the number of regression coefficients
stdlThetaE = sqrt( var(X_E(useindx,:)*betareg(1:length(X_E(1,:)))) ...
                   + mean( (betareg(endbetas+1) + fem(useindx)*betareg(endbetas+2) + blk(useindx)*betareg(endbetas+3) + hsp(useindx)*betareg(endbetas+4)).^2 ) ) ;  %%%%This is the average standard deviation of X_E*betaE+eta_E
stdlThetaL = sqrt( var( X_L(useindx,:)*betareg(length(X_E(1,:))+1:length(X_L(1,:))+length(X_E(1,:))) ) ...
                   + mean( (betareg(endbetas+5) + fem(useindx)*betareg(endbetas+6) + blk(useindx)*betareg(endbetas+7) + hsp(useindx)*betareg(endbetas+8)).^2 ) ) ;  %%%%This is the average standard deviation of X_L*betaL+eta_L
varlThetaEresid = mean( (betareg(endbetas+1) + fem(useindx)*betareg(endbetas+2) + blk(useindx)*betareg(endbetas+3) + hsp(useindx)*betareg(endbetas+4)).^2 ) ;  %%%%This is the average variance of eta_E
varlThetaLresid = mean( (betareg(endbetas+5) + fem(useindx)*betareg(endbetas+6) + blk(useindx)*betareg(endbetas+7) + hsp(useindx)*betareg(endbetas+8)).^2 ) ;  %%%%This is the average variance of eta_L
stdX_E     = std(X_E,'omitnan') ;  %%%%This is a row vector with the standard deviations of the regressors X_E

beta_hatE     = [] ;  %%%%This will store the estimated coefficient values
StDev_EffectE = [] ;  %%%%This will store the standard deviation effect estimates
StdErr_E      = [] ;  %%%%This will store the standard error estimates
cmdnames = 'VarnameE={' ; 
beta_hatE     = [beta_hatE; betareg(1); betareg(2)] ; 
StdErr_E      = [StdErr_E; StdErr(1); StdErr(2)] ; 
StDev_EffectE = [StDev_EffectE; betareg(1)*stdX_E(1)/stdlThetaE ; 
                                betareg(2)*stdX_E(2)/stdlThetaE] ; 
cmdnames = strcat(cmdnames,'''MeanNbhdIncome:'';''UninsuredMinorRate:'';') ; 

if Model==1
    %%%%This is the basic version of the regression model: 
    %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Grade-5, 7 CONST
    beta_hatE     = [beta_hatE; betareg(NbhdVars+1:NbhdVars+5)] ; 
    StdErr_E      = [StdErr_E; StdErr(NbhdVars+1:NbhdVars+5)] ; 
    StDev_EffectE = [StDev_EffectE;
                     (1/stdlThetaE)*betareg(NbhdVars+1:NbhdVars+4); NaN] ; 
    cmdnames     = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';''CONST'';') ; 
elseif Model==2 %%%%This is the basic version of the regression model: 
                %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
    beta_hatE     = [beta_hatE; betareg(NbhdVars+1:NbhdVars+3); ... %fem, black, hispanic
                                betareg(NbhdVars+6); ... %grade-5
                                betareg(NbhdVars+4:NbhdVars+5); ... %Dist2, Dist3
                                betareg(NbhdVars+7)] ; %CONST
                %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
    StdErr_E      = [StdErr_E; StdErr(NbhdVars+1:NbhdVars+3); ... %fem, black, hispanic
                               StdErr(NbhdVars+6); ... %grade-5
                               StdErr(NbhdVars+4:NbhdVars+5); ... %Dist2, Dist3
                               StdErr(NbhdVars+7)] ; %CONST
    StDev_EffectE = [StDev_EffectE; (1/stdlThetaE)*betareg(NbhdVars+1:NbhdVars+3); ... %fem, black, hispanic
                                    (1/stdlThetaE)*betareg(NbhdVars+6); ... %grade-5
                                    (1/stdlThetaE)*betareg(NbhdVars+4:NbhdVars+5); NaN] ;  %Dist2, Dist3, CONST
    cmdnames     = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';',...
                                   '''Dist2'';''Dist3'';''CONST'';') ; 
elseif Model>2 && Model<5 %%%%This is the first 11 variables of regression Models 3 through 4.085: 
                         %%%%theta_E Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                         %%%%                   5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5, 11 CONST
    beta_hatE     = [beta_hatE; betareg(NbhdVars+3:NbhdVars+5); ... %fem, black, hispanic
                                betareg(NbhdVars+8); ... %grade-5
                                betareg(NbhdVars+6:NbhdVars+7); ... %Dist2, Dist3
                                betareg(NbhdVars+1:NbhdVars+2); ...  %helper_adult, helper_peer
                                betareg(NbhdVars+9)] ;  %CONST %%%CONTAINS COEFFICIENT ESTIMATES
    StdErr_E      = [StdErr_E; StdErr(NbhdVars+3:NbhdVars+5); ... %fem, black, hispanic
                               StdErr(NbhdVars+8); ... %grade-5
                               StdErr(NbhdVars+6:NbhdVars+7); ... %Dist2, Dist3
                               StdErr(NbhdVars+1:NbhdVars+2); ...  %helper_adult, helper_peer
                               StdErr(NbhdVars+9)] ;  %CONST %%%CONTAINS COEFFICIENT ESTIMATES
    StDev_EffectE = [StDev_EffectE; (1/stdlThetaE)*betareg(NbhdVars+3:NbhdVars+5); ... %fem, black, hispanic
                                    (1/stdlThetaE)*betareg(NbhdVars+8); ... %grade-5
                                    (1/stdlThetaE)*betareg(NbhdVars+6:NbhdVars+7); ... %Dist2, Dist3
                                    (1/stdlThetaE)*betareg(NbhdVars+1:NbhdVars+2).*stdX_E(NbhdVars+1:NbhdVars+2)'; ...  %helper_adult, helper_peer
                                    NaN] ;  %CONST  %%%%CONTAINS STANDARD DEVIATION EFFECT ESTIMATES
    cmdnames     = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';',...
                                   '''Dist2'';''Dist3'';''#HELPER_ADULT'';''#HELPER_PEER'';''CONST'';') ;
    if Model==4.05 %%%%This is the basic version of the regression model: 
                   %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                   %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
                   %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
        beta_hatE = [beta_hatE; betareg(NbhdVars+10); ... %no_internet
                     betareg(NbhdVars+11:NbhdVars+12)] ;  %mobilefrac, tabletfrac
        StdErr_E  = [StdErr_E; StdErr(NbhdVars+10); ... %no_internet
                     StdErr(NbhdVars+11:NbhdVars+12)] ;  %mobilefrac, tabletfrac
        StDev_EffectE = [StDev_EffectE; (1/stdlThetaE)*betareg(NbhdVars+10); ... %no_internet
                                        (1/stdlThetaE)*betareg(NbhdVars+11:NbhdVars+12).*stdX_E(NbhdVars+11:NbhdVars+12)'] ;  %mobilefrac, tabletfrac
       cmdnames = strcat(cmdnames,'''NO_INTERNET'';''MOBILEFRAC'';''TABLETFRAC'';') ;  

    elseif Model==4.085%%%%This is the basic version of the regression model: 
    %%%%theta_E Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
    %%%%                    5 fem, 6 Black, 7 Hispanic, 8 Dist2, 9 Dist3, 10 Grade-5,
    %%%%                    11 CONST, 12 NO_INTERNET, 13 MOBILEFRAC, 14 TABLETFRAC
    %%%%                    15 InvolvedPar, 16 bigfam, 17 middlechild, 18 youngestchild
        beta_hatE = [beta_hatE; betareg(NbhdVars+10); ... %no_internet
                                betareg(NbhdVars+11:NbhdVars+16)] ;  %MOBILEFRAC, TABLETFRAC, involvedparent, bigfam, middlechild, youngestchild
        StdErr_E = [StdErr_E; StdErr(NbhdVars+10); ... %no_internet
                              StdErr(NbhdVars+11:NbhdVars+16)] ;  %MOBILEFRAC, TABLETFRAC, involvedparent, bigfam, middlechild, youngestchild
        StDev_EffectE = [StDev_EffectE; (1/stdlThetaE)*betareg(NbhdVars+10); ... %no_internet
                                        (1/stdlThetaE)*betareg(NbhdVars+11:NbhdVars+12).*stdX_E(NbhdVars+11:NbhdVars+12)';... %MOBILEFRAC, TABLETFRAC
                                        (1/stdlThetaE)*betareg(NbhdVars+13:NbhdVars+16)] ; %involvedparent, bigfam, middlechild, youngestchild
       cmdnames = strcat(cmdnames,'''NO_INTERNET'';''MOBILEFRAC'';''TABLETFRAC'';''InvolvedPar'';''BigFam'';''MiddleChild'';''YoungestChild'';') ;  
    end
end
    
%%%%NOW WE TACK ON THE VARIANCE PARAMETERS FOR PRODUCTIVITY, ALONG WITH THE PSEUDO-R^2 AND THE SAMPLE SIZE
beta_hatE     = [beta_hatE; 
    ( betareg(endbetas+1) + betareg(endbetas+3)*mean(AchievementData.black==1&AchievementData.fem==0) + betareg(endbetas+4)*mean(AchievementData.hisp==1&AchievementData.fem==0) ); 
    ( betareg(endbetas+1) + betareg(endbetas+2) + betareg(endbetas+3)*mean(AchievementData.black==1&AchievementData.fem==1) + betareg(endbetas+4)*mean(AchievementData.hisp==1&AchievementData.fem==1) ); 
    ( betareg(endbetas+1) + betareg(endbetas+2)*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1)  ); 
    ( betareg(endbetas+1) + betareg(endbetas+3) + betareg(endbetas+2)*mean(AchievementData.black==1&AchievementData.fem==1)  ) ; 
    ( betareg(endbetas+1) + betareg(endbetas+4) + betareg(endbetas+2)*mean(AchievementData.hisp==1&AchievementData.fem==1)  ) ; 
    1-varlThetaEresid/(stdlThetaE^2) ; 
    length(useindx)] ; 
StdErr_E = [StdErr_E;...
            sqrt(StdErr(endbetas+1)^2 + (StdErr(endbetas+3)*mean(AchievementData.black==1&AchievementData.fem==0))^2 + (StdErr(endbetas+4)*mean(AchievementData.hisp==1&AchievementData.fem==0))^2 ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==0)*Vhat(endbetas+1,endbetas+3) + 2*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+1,endbetas+4) ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==0)*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+3,endbetas+4));... %%%%This is the standard error of the average male standard deviation
            sqrt(StdErr(endbetas+1)^2 + StdErr(endbetas+2)^2 ...
                                      + (StdErr(endbetas+3)*mean(AchievementData.black==1&AchievementData.fem==1))^2 + (StdErr(endbetas+4)*mean(AchievementData.hisp==1&AchievementData.fem==1))^2 ...
                                      + 2*Vhat(endbetas+1,endbetas+2) ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==0)*Vhat(endbetas+1,endbetas+3) + 2*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+1,endbetas+4) ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+2,endbetas+3) + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+2,endbetas+4) ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==0)*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+3,endbetas+4));... %%%%This is the standard error of the average female standard deviation
            sqrt(StdErr(endbetas+1)^2 + (StdErr(endbetas+2)*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1))^2 ...
                                      + 2*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1)*Vhat(endbetas+1,endbetas+2));... %%%%This is the standard error of the average White/Asian standard deviation
            sqrt(StdErr(endbetas+1)^2 + StdErr(endbetas+3)^2 + (StdErr(endbetas+2)*mean(AchievementData.black==1&AchievementData.fem==1))^2 ...
                                      + 2*Vhat(endbetas+1,endbetas+3) + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+1,endbetas+2) ...
                                      + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+2,endbetas+3));... %%%%This is the standard error of the average black standard deviation
            sqrt(StdErr(endbetas+1)^2 + StdErr(endbetas+4)^2 + (StdErr(endbetas+2)*mean(AchievementData.hisp==1&AchievementData.fem==1))^2 ...
                                      + 2*Vhat(endbetas+1,endbetas+4) + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+1,endbetas+2) ...
                                      + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+2,endbetas+4));... %%%%This is the standard error of the average hispanic standard deviation
            nan(2,1)] ;
StDev_EffectE = [StDev_EffectE; nan(7,1)] ;
cmdnames     = strcat(cmdnames,'''SigmaE|male'';''SigmaE|female'';''SigmaE|WAO'';''SigmaE|black'';''SigmaE|hispanic'';''R^2 (log(ThetaE))'';''N Obs''};') ;
eval(cmdnames) ;  %%%%THIS CREATES THE VARIABLE VarnamesE, WHICH CONTAINS THE ORDERED LABELS FOR ALL ELEMENTS OF beta_hatE

stdX_L        = std(X_L,'omitnan') ;
beta_hatL     = [] ;
StdErr_L      = [] ;
StDev_EffectL = [] ;
cmdnames = 'VarnameL={' ;
endbetaE = length(X_E(1,:)) ;


beta_hatL = [beta_hatL; betareg(endbetaE+1); betareg(endbetaE+2)] ;
StdErr_L  = [StdErr_L; StdErr(endbetaE+1); StdErr(endbetaE+2)] ;
StDev_EffectL = [StDev_EffectL; betareg(endbetaE+1)*stdX_L(1)/stdlThetaL;
    betareg(endbetaE+2)*stdX_L(2)/stdlThetaL] ;
cmdnames = strcat(cmdnames,'''MeanNbhdIncome:'';''UninsuredMinorRate:'';') ;

if Model==1    %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Grade-5, 7 CONST
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+1:endbetaE+NbhdVars+5)] ; %fem, black, hispanic, grade-5, CONST
    StdErr_L      = [StdErr_L; StdErr(endbetaE+NbhdVars+1:endbetaE+NbhdVars+5)] ; %fem, black, hispanic, grade-5, CONST
    StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+1)/stdlThetaL; ... %fem
        betareg(endbetaE+NbhdVars+2)/stdlThetaL; ... %black
        betareg(endbetaE+NbhdVars+3)/stdlThetaL; ... %hispanic
        betareg(endbetaE+NbhdVars+4)/stdlThetaL; NaN ] ; %grade-5, const
    cmdnames = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';''CONST'';') ;
elseif Model==2  %%%%Variables: 1 NBHD-Inc, 2 NBHD UninsMinors, 3 fem, 4 Black, 5 Hispanic, 6 Dist2, 7 Dist3, 8 Grade-5, 9 CONST
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+1:endbetaE+NbhdVars+3); ... %fem, black, hispanic
                                betareg(endbetaE+NbhdVars+6); ... %grade-5
                                betareg(endbetaE+NbhdVars+4:endbetaE+NbhdVars+5); ... %Dist2, Dist3
                                betareg(endbetaE+NbhdVars+7)] ; %const
    StdErr_L      = [StdErr_L;  StdErr(endbetaE+NbhdVars+1:endbetaE+NbhdVars+3); ... %fem, black, hispanic
                                StdErr(endbetaE+NbhdVars+6); ... %grade-5
                                StdErr(endbetaE+NbhdVars+4:endbetaE+NbhdVars+5); ... %Dist2, Dist3
                                StdErr(endbetaE+NbhdVars+7)] ; %const
    StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+1)/stdlThetaL;... %fem
                                    betareg(endbetaE+NbhdVars+2)/stdlThetaL;... %black
                                    betareg(endbetaE+NbhdVars+3)/stdlThetaL;... %hispanic
                                    betareg(endbetaE+NbhdVars+6)/stdlThetaL;... %grade-5
                                    betareg(endbetaE+NbhdVars+4)/stdlThetaL;... %Dist2
                                    betareg(endbetaE+NbhdVars+5)/stdlThetaL; NaN ] ; %Dist3, const
    cmdnames = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';',...
                               '''DISTRICT_2'';''DISTRICT_3'';''CONST'';') ;
elseif Model==3  %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 HELPER_ADULT, 4 HELPER_PEER, 
                 %%%%                    5 MathFav, 6 MathLeast, 7 EXTRINSIC, 8 INTRINSIC, 9 fem, 
                 %%%%                    10 Black, 11 Hispanic, 12 Dist2, 13 Dist3, 14 Grade-5, 15 CONST
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+7:endbetaE+NbhdVars+9); ... %fem, black, hispanic
                                betareg(endbetaE+NbhdVars+12); ... %grade-5
                                betareg(endbetaE+NbhdVars+10:endbetaE+NbhdVars+11); ... %Dist2, Dist3
                                betareg(endbetaE+NbhdVars+1:endbetaE+NbhdVars+2); ... %helper_adult, helper_peer
                                betareg(endbetaE+NbhdVars+3:endbetaE+NbhdVars+6); ... %mathfav, mathleast, extrinsic, intrinsic
                                betareg(endbetaE+NbhdVars+13)] ; %const
    StdErr_L      = [StdErr_L;  StdErr(endbetaE+NbhdVars+7:endbetaE+NbhdVars+9); ... %fem, black, hispanic
                                StdErr(endbetaE+NbhdVars+12); ... %grade-5
                                StdErr(endbetaE+NbhdVars+10:endbetaE+NbhdVars+11); ... %Dist2, Dist3
                                StdErr(endbetaE+NbhdVars+1:endbetaE+NbhdVars+2); ... %helper_adult, helper_peer
                                StdErr(endbetaE+NbhdVars+3:endbetaE+NbhdVars+6); ... %mathfav, mathleast, extrinsic, intrinsic
                                StdErr(endbetaE+NbhdVars+13)] ; %const
    StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+7)/stdlThetaL;... %fem
                                    betareg(endbetaE+NbhdVars+8)/stdlThetaL;... %black
                                    betareg(endbetaE+NbhdVars+9)/stdlThetaL;... %hispanic
                                    betareg(endbetaE+NbhdVars+12)/stdlThetaL;... %grade-5
                                    betareg(endbetaE+NbhdVars+10)/stdlThetaL;... %Dist2
                                    betareg(endbetaE+NbhdVars+11)/stdlThetaL;... %Dist3
                                    betareg(endbetaE+NbhdVars+1)*stdX_L(NbhdVars+1)/stdlThetaL;... %helper_adult
                                    betareg(endbetaE+NbhdVars+2)*stdX_L(NbhdVars+2)/stdlThetaL;... %helper_peer
                                    betareg(endbetaE+NbhdVars+3)/stdlThetaL;... %mathfav
                                    betareg(endbetaE+NbhdVars+4)/stdlThetaL;... %mathleast
                                    betareg(endbetaE+NbhdVars+5)/stdlThetaL;... %extrinsic
                                    betareg(endbetaE+NbhdVars+6)/stdlThetaL; NaN ] ; %intrinsic, const
    cmdnames = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';',...
                               '''DISTRICT_2'';''DISTRICT_3'';''#HELPER_ADULT'';''#HELPER_PEER'';',...
                               '''MathFav'';''MathLeast'';''EXTRINSIC'';''INTRINSIC'';''CONST'';') ;
elseif Model==4 || Model==4.05 || Model==4.085
                   %%%%HERE ARE THE FIRST 22 VARIABLES FOR EACH OF THESE MODELS:
                   %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
                   %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
                   %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                   %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
                   %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+14:endbetaE+NbhdVars+16); ... %fem, black, hispanic
                                betareg(endbetaE+NbhdVars+19); ... %grade-5
                                betareg(endbetaE+NbhdVars+17:endbetaE+NbhdVars+18); ... %Dist2, Dist3
                                betareg(endbetaE+NbhdVars+4:endbetaE+NbhdVars+7); ... %helper_adult, helper_peer, mathfav, mathleast
                                betareg(endbetaE+NbhdVars+12:endbetaE+NbhdVars+13)] ; %extrinsic, intrinsic
    StdErr_L      = [StdErr_L;  StdErr(endbetaE+NbhdVars+14:endbetaE+NbhdVars+16); ... %fem, black, hispanic
                                StdErr(endbetaE+NbhdVars+19); ... %grade-5
                                StdErr(endbetaE+NbhdVars+17:endbetaE+NbhdVars+18); ... %Dist2, Dist3
                                StdErr(endbetaE+NbhdVars+4:endbetaE+NbhdVars+7); ... %helper_adult, helper_peer, mathfav, mathleast
                                StdErr(endbetaE+NbhdVars+12:endbetaE+NbhdVars+13)] ; %extrinsic, intrinsic
    StDev_EffectL = [StDev_EffectL;... 
        betareg(endbetaE+NbhdVars+14)/stdlThetaL;... %fem
        betareg(endbetaE+NbhdVars+15)/stdlThetaL;... %black
        betareg(endbetaE+NbhdVars+16)/stdlThetaL;... %hispanic
        betareg(endbetaE+NbhdVars+19)/stdlThetaL;... %grade-5
        betareg(endbetaE+NbhdVars+17)/stdlThetaL;... %Dist2
        betareg(endbetaE+NbhdVars+18)/stdlThetaL;... %Dist3
        betareg(endbetaE+NbhdVars+4)*stdX_L(NbhdVars+4)/stdlThetaL;... %helper_adult
        betareg(endbetaE+NbhdVars+5)*stdX_L(NbhdVars+5)/stdlThetaL;... %helper_peer
        betareg(endbetaE+NbhdVars+6)/stdlThetaL;... %mathfav
        betareg(endbetaE+NbhdVars+7)/stdlThetaL;... %mathleast
        betareg(endbetaE+NbhdVars+12)*stdX_L(NbhdVars+12)/stdlThetaL;... %extrinsic
        betareg(endbetaE+NbhdVars+13)*stdX_L(NbhdVars+13)/stdlThetaL] ; %intrinsic
    cmdnames = strcat(cmdnames,'''FEM'';''BLACK'';''HISPANIC'';''Grade-5'';',...
                               '''DISTRICT_2'';''DISTRICT_3'';''#HELPER_ADULT'';''#HELPER_PEER'';',...
                               '''MathFav'';''MathLeast'';''EXTRINSIC'';''INTRINSIC'';') ;
    
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+1);... % #GamingSys 
                                betareg(endbetaE+NbhdVars+2);... % GamingWkDay
                                betareg(endbetaE+NbhdVars+3)] ; %WkDayNet
    StdErr_L      = [StdErr_L;  StdErr(endbetaE+NbhdVars+1);... % #GamingSys 
                                StdErr(endbetaE+NbhdVars+2);... % GamingWkDay
                                StdErr(endbetaE+NbhdVars+3)] ; %WkDayNet
    StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+1)*stdX_L(NbhdVars+1)/stdlThetaL;... % #GamingSys 
                                    betareg(endbetaE+NbhdVars+2)/stdlThetaL;... % GamingWkDay
                                    betareg(endbetaE+NbhdVars+3)*stdX_L(NbhdVars+3)/stdlThetaL] ; %WkDayNet
    cmdnames = strcat(cmdnames,'''#GamingSys'';''GamWkdyAllow'';''WkdyNetRec'';') ;
   
    beta_hatL     = [beta_hatL; betareg(endbetaE+NbhdVars+8:endbetaE+NbhdVars+11); ... %FracSocialSuperv, Sports, music, otheractiv
                                betareg(endbetaE+NbhdVars+20)] ; %Const
    StdErr_L      = [StdErr_L;  StdErr(endbetaE+NbhdVars+8:endbetaE+NbhdVars+11); ... %FracSocialSuperv, Sports, music, otheractiv
                                StdErr(endbetaE+NbhdVars+20)] ; %Const
    StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+8)*stdX_L(NbhdVars+8)/stdlThetaL;... %FracSocialSuperv
                                    betareg(endbetaE+NbhdVars+9)/stdlThetaL;... %Sports
                                    betareg(endbetaE+NbhdVars+10)/stdlThetaL;... %music
                                    betareg(endbetaE+NbhdVars+11)/stdlThetaL; NaN ] ; %otheractiv, const
    cmdnames = strcat(cmdnames,'''FracSocialStruct'';''SPORTS'';''MUSIC'';''OtherActiv'';''CONST'';') ;

    if Model==4.05 %%%%This is the basic version of the regression model: 
                   %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
                   %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
                   %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
                   %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
                   %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
                   %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
        beta_hatL = [beta_hatL; betareg(endbetaE+NbhdVars+21:endbetaE+NbhdVars+25)] ; %timestudy, timescreen, no_internet, mobilefrac, tabletfrac
        StdErr_L  = [StdErr_L; StdErr(endbetaE+NbhdVars+21:endbetaE+NbhdVars+25)] ; %timestudy, timescreen, no_internet, mobilefrac, tabletfrac
        StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+21:endbetaE+NbhdVars+22).*stdX_L(NbhdVars+21:NbhdVars+22)'*(1/stdlThetaL);... %timestudy, timescreen
                                        betareg(endbetaE+NbhdVars+23)*(1/stdlThetaL);... %no_internet
                                        betareg(endbetaE+NbhdVars+24:endbetaE+NbhdVars+25).*stdX_L(NbhdVars+24:NbhdVars+25)'*(1/stdlThetaL)] ; %mobilefrac, tabletfrac
        cmdnames = strcat(cmdnames,'''TIME_STUDY'';''TIME_SCREEN'';''NO_INTERNET'';''MOBILEFRAC'';''TABLETFRAC'';') ;
    elseif Model==4.085 %%%%This is the basic version of the regression model:
    %%%%theta_L Variables:  1 NBHD-Inc, 2 NBHD UninsMinors, 3 #GamingSys, 4 GamingWkDay, 
    %%%%                    5 WkDayNet, 6 HELPER_ADULT, 7 HELPER_PEER, 8 MathFav, 
    %%%%                    9 MathLeast, 10 FracSocialSuperv, 11 Sports, 12 music,
    %%%%                    13 OtherActiv, 14 EXTRINSIC, 15 INTRINSIC, 16 fem, 
    %%%%                    17 Black, 18 Hispanic, 19 Dist2, 20 Dist3, 21 Grade-5, 22 CONST
    %%%%                    23 TimeStudy, 24 TimeScreen, 25 NO_INTERNET, 26 MOBILEFRAC, 27 TABLETFRAC,
    %%%%                    28 InvolvedPar, 29 BigFam, 30 middlechild, 31 youngestchild
        beta_hatL = [beta_hatL; betareg(endbetaE+NbhdVars+21:endbetaE+NbhdVars+29)] ;  %timestudy, timescreen, no_int, mobilefrac, tabletfrac, involvedparent, bigfam, middlechild, youngestchild 
        StdErr_L  = [StdErr_L; StdErr(endbetaE+NbhdVars+21:endbetaE+NbhdVars+29)] ;  %timestudy, timescreen, no_int, mobilefrac, tabletfrac, involvedparent, bigfam, middlechild, youngestchild 
        StDev_EffectL = [StDev_EffectL; betareg(endbetaE+NbhdVars+21:endbetaE+NbhdVars+22).*stdX_L(NbhdVars+21:NbhdVars+22)'*(1/stdlThetaL) ;... %timestudy, timescreen 
                                        betareg(endbetaE+NbhdVars+23)*(1/stdlThetaL) ;... %no_internet 
                                        betareg(endbetaE+NbhdVars+24:endbetaE+NbhdVars+25).*stdX_L(NbhdVars+24:NbhdVars+25)'*(1/stdlThetaL);... %mobilefrac, tabletfrac 
                                        betareg(endbetaE+NbhdVars+26:endbetaE+NbhdVars+29)*(1/stdlThetaL)] ; %involvedparent, bigfam, middlechild, youngestchild 
        cmdnames = strcat(cmdnames,'''TIME_STUDY'';''TIME_SCREEN'';''NO_INTERNET'';''MOBILEFRAC'';''TABLETFRAC'';',...
                                   '''InvolvedPar'';''BigFam'';''MiddleChild'';''YoungestChild'';') ; 
        
    end
end

beta_hatL     = [beta_hatL; 
    ( betareg(endbetas+5) + betareg(endbetas+7)*mean(AchievementData.black==1&AchievementData.fem==0) + betareg(endbetas+8)*mean(AchievementData.hisp==1&AchievementData.fem==0) );
    ( betareg(endbetas+5) + betareg(endbetas+6) + betareg(endbetas+7)*mean(AchievementData.black==1&AchievementData.fem==1) + betareg(endbetas+8)*mean(AchievementData.hisp==1&AchievementData.fem==1) );
    ( betareg(endbetas+5) + betareg(endbetas+6)*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1)  );
    ( betareg(endbetas+5) + betareg(endbetas+7) + betareg(endbetas+6)*mean(AchievementData.black==1&AchievementData.fem==1)  );
    ( betareg(endbetas+5) + betareg(endbetas+8) + betareg(endbetas+6)*mean(AchievementData.hisp==1&AchievementData.fem==1)  ) ;
    1-varlThetaLresid/(stdlThetaL^2);
    length(useindx)] ;
StdErr_L = [StdErr_L; ...
            sqrt( StdErr(endbetas+5)^2 + (StdErr(endbetas+7)*mean(AchievementData.black==1&AchievementData.fem==0))^2 + (betareg(endbetas+8)*mean(AchievementData.hisp==1&AchievementData.fem==0))^2 ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==0)*Vhat(endbetas+5,endbetas+7) ...
                                       + 2*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+5,endbetas+8) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==0)*mean(AchievementData.hisp==1&AchievementData.fem==0)*Vhat(endbetas+7,endbetas+8) ); ... %%%%This is the standard error of the average male standard deviation
            sqrt( StdErr(endbetas+5)^2 + StdErr(endbetas+6)^2 + (StdErr(endbetas+7)*mean(AchievementData.black==1&AchievementData.fem==1))^2 + (betareg(endbetas+8)*mean(AchievementData.hisp==1&AchievementData.fem==1))^2 ...
                                       + 2*Vhat(endbetas+5,endbetas+6) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+5,endbetas+7) ...
                                       + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+5,endbetas+8) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+6,endbetas+7) ...
                                       + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+6,endbetas+8) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==1)*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+7,endbetas+8) ); ... %%%%This is the standard error of the average female standard deviation
            sqrt( StdErr(endbetas+5)^2 + (StdErr(endbetas+6)*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1))^2 ...
                                       + 2*mean(AchievementData.black==0&AchievementData.hisp==0&AchievementData.fem==1)*Vhat(endbetas+5,endbetas+6) ); ... %%%%This is the standard error of the average White/Asian standard deviation
            sqrt( StdErr(endbetas+5)^2 + StdErr(endbetas+7)^2 + (StdErr(endbetas+6)*mean(AchievementData.black==1&AchievementData.fem==1))^2 ...
                                       + 2*Vhat(endbetas+5,endbetas+7) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+5,endbetas+6) ...
                                       + 2*mean(AchievementData.black==1&AchievementData.fem==1)*Vhat(endbetas+6,endbetas+7) ); ... %%%%This is the standard error of the average black standard deviation
            sqrt( StdErr(endbetas+5)^2 + StdErr(endbetas+8)^2 + (StdErr(endbetas+6)*mean(AchievementData.hisp==1&AchievementData.fem==1))^2 ...
                                       + 2*Vhat(endbetas+5,endbetas+8) ...
                                       + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+5,endbetas+6) ...
                                       + 2*mean(AchievementData.hisp==1&AchievementData.fem==1)*Vhat(endbetas+6,endbetas+8) ); ... %%%%This is the standard error of the average black standard deviation
            nan(2,1)] ;
StDev_EffectL = [StDev_EffectL; nan(7,1)] ;
cmdnames     = strcat(cmdnames,'''SigmaL|male'';''SigmaL|female'';''SigmaL|WAO'';''SigmaL|black'';''SigmaL|hispanic'';''R^2 (log(ThetaL))'';''N Obs''};') ;
eval(cmdnames) ;

format long


disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%                 TABLES 2 & 3:             %%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')

betaE = table(VarnameE,beta_hatE,StDev_EffectE,StdErr_E) ; 
disp(betaE(:,:))


betaL = table(VarnameL,beta_hatL,StDev_EffectL,StdErr_L) ; 
disp(betaL(:,:))






%%

disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%  STAGE 4 ESTIMATION:')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')
disp('-')

K_E        = length(X_E(1,:)) ;
K_L        = length(X_L(1,:)) ;
Q          = AchievementData.Q ; Q = Q(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
T          = AchievementData.T ; T = T(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
D2         = AchievementData.Dist2 ; D2 = D2(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
D3         = AchievementData.Dist3 ; D3 = D3(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
grade5     = AchievementData.grade5 ; grade5 = grade5(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
SSEInc     = X_E(:,1) ; SSEInc = SSEInc(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
SSEIns     = X_E(:,2) ; SSEIns = SSEIns(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
con        = con(useindx) ;
wts        = wts(useindx) ;
fem        = fem(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
blk        = blk(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
hsp        = hsp(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
Y_E        = Y_E(useindx) ;
Y_L        = Y_L(useindx) ;
X_E        = X_E(useindx,:) ;
X_L        = X_L(useindx,:) ;
%%%%The following several lines compute EB forecasts of pre-test scores
%%%%using the IRT setup from Appendix A.5.1:
precorrect     = (AchievementData.pretest + 9 -0.25*AchievementData.preblank)/1.25 ; 
kappapre       = precorrect/36 ;  kappapre(kappapre==0) = 1/(36*2) ; kappapre(kappapre==1) = 1-1/(36*2) ;
sig_kappapre   = kappapre.*(1-kappapre)/36 ; 
SIG_kappapre   = mean(sig_kappapre,'omitnan') ;
lambdakappapre = SIG_kappapre./(SIG_kappapre + sig_kappapre) ; 
MU_kappapre    = mean(kappapre,'omitnan') ; 
kappapreEB     = lambdakappapre.*kappapre + (1-lambdakappapre)*MU_kappapre ;

pre             = kappapreEB ; pre = pre(useindx) ;
%%%%The following several lines compute EB forecasts of post-test scores
%%%%using the IRT setup from Appendix A.5.1:
postcorrect     = (AchievementData.posttest + 9 -0.25*AchievementData.postblank)/1.25 ; 
kappapost       = postcorrect/36 ;  kappapost(kappapost==0) = 1/(36*2) ; kappapost(kappapost==1) = 1-1/(36*2) ;
sig_kappapost   = kappapost.*(1-kappapost)/36 ; 
SIG_kappapost   = mean(sig_kappapost,'omitnan') ;
lambdakappapost = SIG_kappapost./(SIG_kappapost + sig_kappapost) ; 
MU_kappapost    = mean(kappapost,'omitnan') ; 
kappapostEB     = lambdakappapost.*kappapost + (1-lambdakappapost)*MU_kappapost ;

post            = kappapostEB ; post = post(useindx) ; 
Deltakappa      = kappapost - kappapre ;  Deltakappa = Deltakappa(useindx) ; 
DeltakappaEB    = kappapostEB - kappapreEB ; DeltakappaEB = DeltakappaEB(useindx) ;  

disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%       TABLE 5 DESCRIPTIVE STATISTICS      %%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
disp('--------------------- PRE-TEST --------------------------------------')
meanSD_kappapre    = [mean(kappapre(~isnan(kappapost))*36) std(kappapre(~isnan(kappapost))*36) std(kappapreEB(~isnan(kappapost))*36) sum(~isnan(kappapost))] %#ok
meanSD_kappapreFEM = [mean(kappapre(~isnan(kappapost)&g2.female==1)*36) std(kappapre(~isnan(kappapost)&g2.female==1)*36) std(kappapreEB(~isnan(kappapost)&g2.female==1)*36) sum(~isnan(kappapost)&(g2.female==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapreMAL = [mean(kappapre(~isnan(kappapost)&g2.female==0)*36) std(kappapre(~isnan(kappapost)&g2.female==0)*36) std(kappapreEB(~isnan(kappapost)&g2.female==0)*36) sum(~isnan(kappapost)&(g2.female==0))/sum(~isnan(kappapost))] %#ok
meanSD_kappapreBLK = [mean(kappapre(~isnan(kappapost)&g2.black==1)*36) std(kappapre(~isnan(kappapost)&g2.black==1)*36) std(kappapreEB(~isnan(kappapost)&g2.black==1)*36) sum(~isnan(kappapost)&(g2.black==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapreHSP = [mean(kappapre(~isnan(kappapost)&g2.hispanic==1)*36) std(kappapre(~isnan(kappapost)&g2.hispanic==1)*36) std(kappapreEB(~isnan(kappapost)&g2.hispanic==1)*36) sum(~isnan(kappapost)&(g2.hispanic==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapreWAO = [mean(kappapre(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) std(kappapre(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) std(kappapreEB(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) sum(~isnan(kappapost)&(g2.black==0&g2.hispanic==0))/sum(~isnan(kappapost))] %#ok
disp('--------------------- POST-TEST --------------------------------------')
disp('-')
meanSD_kappapost    = [mean(kappapost(~isnan(kappapost))*36) std(kappapost(~isnan(kappapost))*36) std(kappapostEB(~isnan(kappapost))*36) sum(~isnan(kappapost))] %#ok
meanSD_kappapostFEM = [mean(kappapost(~isnan(kappapost)&g2.female==1)*36) std(kappapost(~isnan(kappapost)&g2.female==1)*36) std(kappapostEB(~isnan(kappapost)&g2.female==1)*36) sum(~isnan(kappapost)&(g2.female==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapostMAL = [mean(kappapost(~isnan(kappapost)&g2.female==0)*36) std(kappapost(~isnan(kappapost)&g2.female==0)*36) std(kappapostEB(~isnan(kappapost)&g2.female==0)*36) sum(~isnan(kappapost)&(g2.female==0))/sum(~isnan(kappapost))] %#ok
meanSD_kappapostBLK = [mean(kappapost(~isnan(kappapost)&g2.black==1)*36) std(kappapost(~isnan(kappapost)&g2.black==1)*36) std(kappapostEB(~isnan(kappapost)&g2.black==1)*36) sum(~isnan(kappapost)&(g2.black==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapostHSP = [mean(kappapost(~isnan(kappapost)&g2.hispanic==1)*36) std(kappapost(~isnan(kappapost)&g2.hispanic==1)*36) std(kappapostEB(~isnan(kappapost)&g2.hispanic==1)*36) sum(~isnan(kappapost)&(g2.hispanic==1))/sum(~isnan(kappapost))] %#ok
meanSD_kappapostWAO = [mean(kappapost(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) std(kappapost(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) std(kappapostEB(~isnan(kappapost)&g2.black==0&g2.hispanic==0)*36) sum(~isnan(kappapost)&(g2.black==0&g2.hispanic==0))/sum(~isnan(kappapost))] %#ok
disp('--------------------- TEST SCORE CHANGES  --------------------------------------')
disp('-')
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa))) ;
meanSD_Deltakappa    = [mean(Deltakappa*36,'omitnan') std(Deltakappa*36,'omitnan') std(DeltakappaEB*36,'omitnan') temp sum(~isnan(Deltakappa))] %#ok
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.female==1)) ;
meanSD_DeltakappaFEM = [mean(Deltakappa(g2.female==1)*36,'omitnan') std(Deltakappa(g2.female==1)*36,'omitnan') std(DeltakappaEB(g2.female==1)*36,'omitnan') temp] %#ok
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.female==0)) ;
meanSD_DeltakappaMAL = [mean(Deltakappa(g2.female==0)*36,'omitnan') std(Deltakappa(g2.female==0)*36,'omitnan') std(DeltakappaEB(g2.female==0)*36,'omitnan') temp] %#ok
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.black==1)) ;
meanSD_DeltakappaBLK = [mean(Deltakappa(g2.black==1)*36,'omitnan') std(Deltakappa(g2.black==1)*36,'omitnan') std(DeltakappaEB(g2.black==1)*36,'omitnan') temp] %#ok
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.hispanic==1)) ;
meanSD_DeltakappaHSP = [mean(Deltakappa(g2.hispanic==1)*36,'omitnan') std(Deltakappa(g2.hispanic==1)*36,'omitnan') std(DeltakappaEB(g2.hispanic==1)*36,'omitnan') temp] %#ok
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&(g2.black==0&g2.hispanic==0))) ;
meanSD_DeltakappaWAO = [mean(Deltakappa(g2.black==0&g2.hispanic==0)*36,'omitnan') std(Deltakappa(g2.black==0&g2.hispanic==0)*36,'omitnan') std(DeltakappaEB(g2.black==0&g2.hispanic==0)*36,'omitnan') temp] %#ok


disp('--------------------- TEST SCORE CHANGES BY ACTIVE STATUS  --------------------------------------')
disp('-')
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.total_pass>=2)) ; 
meanSD_DeltakappaACTIVE = [mean(Deltakappa(g2.total_pass>=2)*36,'omitnan') std(Deltakappa(g2.total_pass>=2)*36,'omitnan') std(DeltakappaEB(g2.total_pass>=2)*36,'omitnan') temp] %#ok 
[~,temp] = ttest(Deltakappa(~isnan(Deltakappa)&g2.total_pass<2)) ; 
meanSD_DeltakappaMARGINACT = [mean(Deltakappa(g2.total_pass<2)*36,'omitnan') std(Deltakappa(g2.total_pass<2)*36,'omitnan') std(DeltakappaEB(g2.total_pass<2)*36,'omitnan') temp] %#ok 

[H0_Act_vs_MA_equal,P_Act_vs_MA_equal] = ttest2(Deltakappa(~isnan(Deltakappa)&g2.total_pass>=2),Deltakappa(~isnan(Deltakappa)&g2.total_pass<2)) %#ok
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')





preblank   = AchievementData.preblank ; preblank = preblank(useindx) ;
postblank  = AchievementData.postblank ; postblank = postblank(useindx) ;
Deltablank = (36-postblank) - (36-preblank) ;
BETA_E     = betareg(1:K_E) ;
BETA_L     = betareg(K_E+1:K_E+K_L) ;
muE        = 0 ;
muL        = 0 ;
sigmaE     = betareg(K_E+K_L+1) ;
sigmaE_f   = betareg(K_E+K_L+2) ;
sigmaE_blk = betareg(K_E+K_L+3) ;
sigmaE_hsp = betareg(K_E+K_L+4) ;
SIGMAe     = sigmaE + sigmaE_f*fem + sigmaE_blk*blk + sigmaE_hsp*hsp ;
sigmaL     = betareg(K_E+K_L+5) ;
sigmaL_f   = betareg(K_E+K_L+6) ;
sigmaL_blk = betareg(K_E+K_L+7) ;
sigmaL_hsp = betareg(K_E+K_L+8) ;
SIGMAl     = sigmaL + sigmaL_f*fem + sigmaL_blk*blk + sigmaL_hsp*hsp ;
rhoEL      = betareg(K_E+K_L+9) ;
rhoEL_f    = betareg(K_E+K_L+10) ;
rhoEL_blk  = betareg(K_E+K_L+11) ;
rhoEL_hsp  = betareg(K_E+K_L+12) ;
RHOel      = rhoEL + rhoEL_f*fem + rhoEL_blk*blk + rhoEL_hsp*hsp ;




LPointEst.ThetaL = nan(size(LPointEst.ThetaLACT(:,1))) ;
for ii=1:length(LPointEst.ThetaL)
    LPointEst.ThetaL(ii) = AchievementData.ThetaL(AchievementData.id==LPointEst.Activeids(ii)) ; 
end





%%%%BSSamp{s,1} contains the student id's (relative to the original dataset) for the re-sampled students in the s-th bootstrap sample.
%%%%tempid = BSSamp{s,1}(ii), ii=1,...,1676 contains the student id that will map the above into the estimated samples of types and covariates.
%%%%indxEL  = find(ThetaEThetaLBSid==tempid) is the row of ThetaEBS and ThetaLBS for student tempid
%%%%ThetaEBS(indxEL,s) contains the ThetaE type for student id tempid in bootstrap sample s
%%%%ThetaLBS(indxEL,s) contains the ThetaE type for student id tempid in bootstrap sample s
%%%%indxX   = find(Xid==tempid) is the row of X_E and X_L for student tempid
%%%%X_E(indxX,:) contains the log(ThetaE)-covariates for student id tempid
%%%%X_L(indxX,:) contains the log(ThetaL)-covariates for student id tempid
%%%%Q(indxX) contains the total output level for student id tempid 
%%%%con(indxX) contains the contract for student id tempid 
%%%%fem(indxX) contains the gender for student id tempid 
%%%%blk(indxX) contains the black binary for student id tempid 
%%%%hsp(indxX) contains the hispanic binary for student id tempid 
%%%%wts(indxX) contains the sample weights for student id tempid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LEARNING PRODUCTION FUNCTION REGRESSION: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%Step 1: Compute E[theta|X,Q<2] for each person in the selected-out set. 
%%%%NOTE THAT IN ORDER TO DO SO, WE NEED TO INTEGRATE OVER THE CONDITIONAL JOINT
%%%%DISTRIBUTION OF (etaE,etaL), GIVEN THAT (thetaE,thetaL) WAS SUCH THAT Q<2.  
%%%%WE DO THIS INTEGRATION VIA MONTE CARLO METHODS.
tic 
rng(round(pi*1000000)) ; 
NsimZ = 1000000 ; 
SimZ  = [randn(NsimZ,1) randn(NsimZ,1)] ; 
Y_Eimpute = Y_E ; 
Y_Limpute = Y_L ; 
parfor ii=1:length(Q) 
    if Q(ii)<2 
        M = [muE muL] ;
        S = [SIGMAe(ii)^2, RHOel(ii)*SIGMAe(ii)*SIGMAl(ii); ...
             RHOel(ii)*SIGMAe(ii)*SIGMAl(ii), SIGMAl(ii)^2] ; 
        
        %%%%The following code transforms a 2-D vector of iid standard normal values into
        %%%%general random variables with mean vector M and variance-covariance matrix S.
        %%%%STEP 1: compute the Cholesky factorization of the var-cov matrix:
        V = chol(S) ;  %%%%NOTE: S = V'*V
        %%%%STEP 2: here is how we convert all of the Z-values in SimZ to general Normal RVs
        SimNorm = ( V'*SimZ' + M'*ones(1,NsimZ) )' ;
        %%%%STEP 3: add the regressors to the simulated residuals with the appropriate
        %%%%distribution from above and we have complete simulations from (ThetaE,ThetaL) space.
        SimNorm(:,1) = SimNorm(:,1) + X_E(ii,:)*BETA_E ;
        SimNorm(:,2) = SimNorm(:,2) + X_L(ii,:)*BETA_L ;
        if con(ii)==1
            outflag = ( SimNorm(:,2) > interp1(log(LPointEst.domE),log(LPointEst.threshold1),SimNorm(:,1),'linear','extrap') ) ; %#ok
            Y_Eimpute(ii) = exp(mean(SimNorm(outflag,1))) ;
            Y_Limpute(ii) = exp(mean(SimNorm(outflag,2))) ;
        elseif con(ii)==2
            outflag = ( SimNorm(:,2) > interp1(log(LPointEst.domE),log(LPointEst.threshold2),SimNorm(:,1),'linear','extrap') ) ;
            Y_Eimpute(ii) = exp(mean(SimNorm(outflag,1))) ;
            Y_Limpute(ii) = exp(mean(SimNorm(outflag,2))) ;
        elseif con(ii)==3
            outflag = ( SimNorm(:,2) > interp1(log(LPointEst.domE),log(LPointEst.threshold3),SimNorm(:,1),'linear','extrap') ) ;
            Y_Eimpute(ii) = exp(mean(SimNorm(outflag,1))) ;
            Y_Limpute(ii) = exp(mean(SimNorm(outflag,2))) ;
        end
    end
end
toc

%%%%Here we compute the individual bootstrapped variance of each element of Y_Eimpute and Y_Limpute,
%%%%so that we can use that information to compute an Empirical Bayes shrunk forecast:
BootstrapY_impute = 0 ;
if BootstrapY_impute==1

    %%%%START BY EXPANDING THE TYPE BOOTSTRAP MATRICES TO CONTAIN ROWS FOR THE ENTIRE SAMPLE:
    ThetaEBSoriginal = ThetaEBS ; %%%%The versions of ThetaEBS and ThetaLBS computed above only contain numeric entries for active students.
    ThetaLBSoriginal = ThetaLBS ;
    ThetaEBS     = nan(length(useindx),NBS) ;
    ThetaLBS     = nan(length(useindx),NBS) ;
    useids       = AchievementData.id(useindx) ; %%%%This is the set of student ids from the full samle
    %%%%IMPORTANT NOTE: LPointEst.student_id is the list of all students with Q>=2
    %%%%or  from contracts 1 and 2 with Q==1.  This list of ids determines the
    %%%%identity of the person in the rows of ThetaLBS and ThetaEBS, as computed in
    %%%%ProcessBootstraps.m.
    for ii=1:length(useids)
        if sum(useids(ii)==LPointEst.student_id)>0
            ThetaEBS(ii,:) = ThetaEBSoriginal(LPointEst.student_id==useids(ii),:) ;
            ThetaLBS(ii,:) = ThetaLBSoriginal(LPointEst.student_id==useids(ii),:) ;
        end
    end

    GenParametricBootstrapBetareg = 0 ;
    if GenParametricBootstrapBetareg==1
        run ParametricTobitBootstrap_ThetaELse
    else
        if Model==4.05 
            load('.\STAGE4\betaregPBS_Model4point05.mat') 
        elseif Model==4.085 
            load('.\STAGE4\betaregPBS_Model4point085.mat') 
        end
    end
    Xid              = AchievementData.id ;
    Y_EimputeBS = nan(length(Y_Eimpute),NBS) ;   %%%%These two matrices will contain 1676 rows (one for each student) and NBS columns (one for each bootstrap).
    Y_LimputeBS = nan(length(Y_Limpute),NBS) ;   %%%%The rows corresponding to active students will contain bootstrapped type estimates computed using the first-
                                                 %%%%stage structural estimates.  The rows corresponding to inactive students will contiain bootstrapped type 
                                                 %%%%estimates computed using the Tobit barametric bootstrap samples of betareg (stored in betaregPBS).
    %%%%Here we do some preliminary variable definitions to avoid
    %%%%referencing a broadcast variable in the parfor loop below
    ldomEtemp       = log(LPointEst.domE) ; 
    lthreshold1temp = log(LPointEst.threshold1) ; 
    lthreshold2temp = log(LPointEst.threshold2) ; 
    lthreshold3temp = log(LPointEst.threshold3) ; 
    for s=1:NBS
        tic
        s %#ok
        Y_EimputeBS(:,s) = ThetaEBS(useindx,s) ;  %%%%Note that the bootstrapped type estimates come pre-populated for active students; all other entries are nans.
        Y_LimputeBS(:,s) = ThetaLBS(useindx,s) ;  %%%%Note that the bootstrapped type estimates come pre-populated for active students; all other entries are nans.
        betaEtemp        = betaregPBS(s,1:K_E)' ;
        betaLtemp        = betaregPBS(s,K_E+1:K_E+K_L)' ;
        sigmaEtemp       = betaregPBS(s,K_E+K_L+1) ;
        sigmaE_ftemp     = betaregPBS(s,K_E+K_L+2) ;
        sigmaE_blktemp   = betaregPBS(s,K_E+K_L+3) ;
        sigmaE_hsptemp   = betaregPBS(s,K_E+K_L+4) ;
        SIGMAetemp       = sigmaEtemp + sigmaE_ftemp*fem + sigmaE_blktemp*blk + sigmaE_hsptemp*hsp ;
        sigmaLtemp       = betaregPBS(s,K_E+K_L+5) ;
        sigmaL_ftemp     = betaregPBS(s,K_E+K_L+6) ;
        sigmaL_blktemp   = betaregPBS(s,K_E+K_L+7) ;
        sigmaL_hsptemp   = betaregPBS(s,K_E+K_L+8) ;
        SIGMAltemp       = sigmaLtemp + sigmaL_ftemp*fem + sigmaL_blktemp*blk + sigmaL_hsptemp*hsp ;
        rhoELtemp        = betaregPBS(s,K_E+K_L+9) ;
        rhoEL_ftemp      = betaregPBS(s,K_E+K_L+10) ;
        rhoEL_blktemp    = betaregPBS(s,K_E+K_L+11) ;
        rhoEL_hsptemp    = betaregPBS(s,K_E+K_L+12) ;
        RHOeltemp        = rhoELtemp + rhoEL_ftemp*fem + rhoEL_blktemp*blk + rhoEL_hsptemp*hsp ;
               
        parfor ii=1:length(Q)
            if Q(ii)<2
                M = [muE muL] ;
                S = [SIGMAetemp(ii)^2, RHOeltemp(ii)*SIGMAetemp(ii)*SIGMAltemp(ii); ...
                    RHOeltemp(ii)*SIGMAetemp(ii)*SIGMAltemp(ii), SIGMAltemp(ii)^2] ;
                
                %%%%The following code transforms a 2-D vector of iid standard normal values into
                %%%%general random variables with mean vector M and variance-covariance matrix S.
                %%%%STEP 1: compute the Cholesky factorization of the var-cov matrix:
                V = chol(S) ;  %%%%NOTE: S = V'*V
                %%%%STEP 2: here is how we convert all of the Z-values in SimZ to general Normal RVs
                SimNorm = ( V'*SimZ' + M'*ones(1,NsimZ) )' ;
                %%%%STEP 3: add the regressors to the simulated residuals with the appropriate
                %%%%distribution from above and we have complete simulations from (ThetaE,ThetaL) space.
                SimNorm(:,1) = SimNorm(:,1) + X_E(ii,:)*betaEtemp ;
                SimNorm(:,2) = SimNorm(:,2) + X_L(ii,:)*betaLtemp ;
                if con(ii)==1
                    outflag = ( SimNorm(:,2) > interp1(ldomEtemp,lthreshold1temp,SimNorm(:,1),'linear','extrap') ) ;
                    Y_EimputeBS(ii,s) = exp(mean(SimNorm(outflag,1))) ;
                    Y_LimputeBS(ii,s) = exp(mean(SimNorm(outflag,2))) ;
                elseif con(ii)==2
                    outflag = ( SimNorm(:,2) > interp1(ldomEtemp,lthreshold2temp,SimNorm(:,1),'linear','extrap') ) ;
                    Y_EimputeBS(ii,s) = exp(mean(SimNorm(outflag,1))) ;
                    Y_LimputeBS(ii,s) = exp(mean(SimNorm(outflag,2))) ;
                elseif con(ii)==3
                    outflag = ( SimNorm(:,2) > interp1(ldomEtemp,lthreshold3temp,SimNorm(:,1),'linear','extrap') ) ;
                    Y_EimputeBS(ii,s) = exp(mean(SimNorm(outflag,1))) ;
                    Y_LimputeBS(ii,s) = exp(mean(SimNorm(outflag,2))) ;
                end
            end
        end
        toc
    end
    if Model==4.05
        save('.\STAGE4\BootstrapY_imputeModel4point05.mat','Y_EimputeBS','Y_LimputeBS')
    elseif Model==4.085
        save('.\STAGE4\BootstrapY_imputeModel4point085.mat','Y_EimputeBS','Y_LimputeBS')
    end
    clear betaEtemp betaLtemp ldomEtemp lthreshold1temp lthreshold2temp lthreshold3temp
    clear SIGMAetemp sigmaEtemp sigmaE_ftemp sigmaE_blktemp sigmaE_hsptemp
    clear SIGMAltemp sigmaLtemp sigmaL_ftemp sigmaL_blktemp sigmaL_hsptemp
    clear RHOeltemp rhoELtemp rhoEL_ftemp rhoEL_blktemp rhoEL_hsptemp
else
    if Model==4.05
        load('.\STAGE4\BootstrapY_imputeModel4point05.mat')
    elseif Model==4.085
        load('.\STAGE4\BootstrapY_imputeModel4point085.mat')
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  THIS BLOCK OF CODE COMPUTES MEAN (ACROSS ALL STUDENTS IN THE SAMPLE) VARIANCES AND 
%%%%  CORRELATIONS OF XjBETAj, etaj, AND (thetaP,thetaM), WITH PARAMETRIC BOOTSTRAPPED ESTIMATES OF
%%%%  THE STANDARD ERRORS FOR THESE NUMBERS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic 
rng(round(pi*1000000)) ; 
NsimZ = 10000 ; 
SimZ  = [randn(NsimZ,1) randn(NsimZ,1)] ; 
Cov_eta    = nan(length(Q),1) ; 
Var_eta_E  = nan(length(Q),1) ; 
Var_eta_L  = nan(length(Q),1) ; 
XBeta_E    = X_E*BETA_E ; 
XBeta_L    = X_L*BETA_L ; 
temp       = cov(XBeta_E,XBeta_L) ;
Cov_XBeta  = temp(1,2) ; 
Var_XBeta  = [temp(1,1), temp(2,2)] ; 
for ii=1:length(Q)
    M = [muE muL] ;
    S = [SIGMAe(ii)^2, RHOel(ii)*SIGMAe(ii)*SIGMAl(ii); ...
        RHOel(ii)*SIGMAe(ii)*SIGMAl(ii), SIGMAl(ii)^2] ;

    %%%%The following code transforms a 2-D vector of iid standard normal values into
    %%%%general random variables with mean vector M and variance-covariance matrix S.
    %%%%STEP 1: compute the Cholesky factorization of the var-cov matrix:
    V = chol(S) ;  %%%%NOTE: S = V'*V
    %%%%STEP 2: here is how we convert all of the Z-values in SimZ to general Normal RVs
    eta = ( V'*SimZ' + M'*ones(1,NsimZ) )' ;
    %%%%STEP 3: add the regressors to the simulated residuals with the appropriate
    %%%%distribution from above and we have complete simulations from (ThetaE,ThetaL) space.
    
    temp           = cov(eta) ;
    Cov_eta(ii)    = temp(1,2) ; 
    Var_eta_E(ii)  = temp(1,1) ; 
    Var_eta_L(ii)  = temp(2,2) ;

end
STDEVs_E = [sqrt(Var_XBeta(1)), sqrt(mean(Var_eta_E)), sqrt(mean(Var_eta_E) + Var_XBeta(1))]  ;
STDEVs_L = [sqrt(Var_XBeta(2)), sqrt(mean(Var_eta_L)), sqrt(mean(Var_eta_L) + Var_XBeta(2))]  ;
RHO      = [Cov_XBeta/(STDEVs_E(1)*STDEVs_L(1)), mean(Cov_eta)/(STDEVs_E(2)*STDEVs_L(2)), (mean(Cov_eta) + Cov_XBeta)/(STDEVs_E(3)*STDEVs_L(3))] ;
toc
STDEVs_EPBS1 = nan(400,1) ;
STDEVs_EPBS2 = nan(400,1) ;
STDEVs_EPBS3 = nan(400,1) ;
STDEVs_LPBS1 = nan(400,1) ;
STDEVs_LPBS2 = nan(400,1) ;
STDEVs_LPBS3 = nan(400,1) ;
RHOPBS1      = nan(400,1) ;
RHOPBS2      = nan(400,1) ;
RHOPBS3      = nan(400,1) ;


if Model==4.05
    load('.\STAGE4\betaregPBS_Model4point05.mat')
elseif Model==4.085
    load('.\STAGE4\betaregPBS_Model4point085.mat')
end

parfor ss=1:400
    BETA_EPBS     = betaregPBS(ss,1:K_E)' ;
    BETA_LPBS     = betaregPBS(ss,K_E+1:K_E+K_L)' ;
    muEPBS        = 0 ;
    muLPBS        = 0 ;
    sigmaEPBS     = betaregPBS(ss,K_E+K_L+1) ;
    sigmaE_fPBS   = betaregPBS(ss,K_E+K_L+2) ;
    sigmaE_blkPBS = betaregPBS(ss,K_E+K_L+3) ;
    sigmaE_hspPBS = betaregPBS(ss,K_E+K_L+4) ;
    SIGMAePBS     = sigmaEPBS + sigmaE_fPBS*fem + sigmaE_blkPBS*blk + sigmaE_hspPBS*hsp ;
    sigmaLPBS     = betaregPBS(ss,K_E+K_L+5) ;
    sigmaL_fPBS   = betaregPBS(ss,K_E+K_L+6) ;
    sigmaL_blkPBS = betaregPBS(ss,K_E+K_L+7) ;
    sigmaL_hspPBS = betaregPBS(ss,K_E+K_L+8) ;
    SIGMAlPBS     = sigmaLPBS + sigmaL_fPBS*fem + sigmaL_blkPBS*blk + sigmaL_hspPBS*hsp ;
    rhoELPBS      = betaregPBS(ss,K_E+K_L+9) ;
    rhoEL_fPBS    = betaregPBS(ss,K_E+K_L+10) ;
    rhoEL_blkPBS  = betaregPBS(ss,K_E+K_L+11) ;
    rhoEL_hspPBS  = betaregPBS(ss,K_E+K_L+12) ;
    RHOelPBS      = rhoELPBS + rhoEL_fPBS*fem + rhoEL_blkPBS*blk + rhoEL_hspPBS*hsp ;


    rng(round(pi*1000000)) ;
    NsimZ = 10000 ;
    SimZ  = [randn(NsimZ,1) randn(NsimZ,1)] ;
    Cov_etaPBS    = nan(length(Q),1) ;
    Var_eta_EPBS  = nan(length(Q),1) ;
    Var_eta_LPBS  = nan(length(Q),1) ;
    XBeta_EPBS    = X_E*BETA_EPBS ;
    XBeta_LPBS    = X_L*BETA_LPBS ;
    temp          = cov(XBeta_EPBS,XBeta_LPBS) ;
    Cov_XBetaPBS  = temp(1,2) ;
    Var_XBetaPBS  = [temp(1,1), temp(2,2)] ;
    for ii=1:length(Q)
        M = [muEPBS muLPBS] ;
        S = [SIGMAePBS(ii)^2, RHOelPBS(ii)*SIGMAePBS(ii)*SIGMAlPBS(ii); ...
            RHOelPBS(ii)*SIGMAePBS(ii)*SIGMAlPBS(ii), SIGMAlPBS(ii)^2] ;

        %%%%The following code transforms a 2-D vector of iid standard normal values into
        %%%%general random variables with mean vector M and variance-covariance matrix S.
        %%%%STEP 1: compute the Cholesky factorization of the var-cov matrix:
        V = chol(S) ;  %%%%NOTE: S = V'*V
        %%%%STEP 2: here is how we convert all of the Z-values in SimZ to general Normal RVs
        eta = ( V'*SimZ' + M'*ones(1,NsimZ) )' ;
        %%%%STEP 3: add the regressors to the simulated residuals with the appropriate
        %%%%distribution from above and we have complete simulations from (ThetaE,ThetaL) space.

        temp              = cov(eta) ;
        Cov_etaPBS(ii)    = temp(1,2) ;
        Var_eta_EPBS(ii)  = temp(1,1) ;
        Var_eta_LPBS(ii)  = temp(2,2) ;

    end
    STDEVs_EPBS = [sqrt(Var_XBetaPBS(1)), sqrt(mean(Var_eta_EPBS)), sqrt(mean(Var_eta_EPBS) + Var_XBetaPBS(1))] ;  
    STDEVs_EPBS1(ss) = STDEVs_EPBS(1) ; 
    STDEVs_EPBS2(ss) = STDEVs_EPBS(2) ; 
    STDEVs_EPBS3(ss) = STDEVs_EPBS(3) ; 
    STDEVs_LPBS = [sqrt(Var_XBetaPBS(2)), sqrt(mean(Var_eta_LPBS)), sqrt(mean(Var_eta_LPBS) + Var_XBetaPBS(2))] ;  
    STDEVs_LPBS1(ss) = STDEVs_LPBS(1) ; 
    STDEVs_LPBS2(ss) = STDEVs_LPBS(2) ; 
    STDEVs_LPBS3(ss) = STDEVs_LPBS(3) ; 
    RHOPBS      = [Cov_XBetaPBS/(STDEVs_EPBS(1)*STDEVs_LPBS(1)), mean(Cov_etaPBS)/(STDEVs_EPBS(2)*STDEVs_LPBS(2)), (mean(Cov_etaPBS) + Cov_XBetaPBS)/(STDEVs_EPBS(3)*STDEVs_LPBS(3))] ; 
    RHOPBS1(ss) = RHOPBS(1) ; 
    RHOPBS2(ss) = RHOPBS(2) ; 
    RHOPBS3(ss) = RHOPBS(3) ; 
end
STDEVs_EPBS = sort([STDEVs_EPBS1 STDEVs_EPBS2 STDEVs_EPBS3]) - mean([STDEVs_EPBS1 STDEVs_EPBS2 STDEVs_EPBS3]) + STDEVs_E ; 
STDEVs_LPBS = sort([STDEVs_LPBS1 STDEVs_LPBS2 STDEVs_LPBS3]) - mean([STDEVs_LPBS1 STDEVs_LPBS2 STDEVs_LPBS3]) + STDEVs_L ; 
RHOPBS      = sort([RHOPBS1 RHOPBS2 RHOPBS3]) - mean([RHOPBS1 RHOPBS2 RHOPBS3]) + RHO ; 

disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%       TABLE 4 CORRELATION RESULTS:        %%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')
format shortg
STDResults = {'','SD(XBetaE)','SD(XBetaL)','SD(etaE)','SD(etaL)','SD(THETA_E)','SD(THETA_L)'; ...
              'PE',STDEVs_E(1),STDEVs_L(1),STDEVs_E(2),STDEVs_L(2),STDEVs_E(3),STDEVs_L(3); ...
              'StdErr',std(STDEVs_EPBS(:,1)),std(STDEVs_LPBS(:,1)),std(STDEVs_EPBS(:,2)),std(STDEVs_LPBS(:,2)),std(STDEVs_EPBS(:,3)),std(STDEVs_LPBS(:,3)); ...
              'CI95upper',STDEVs_EPBS(390,1),STDEVs_LPBS(390,1),STDEVs_EPBS(390,2),STDEVs_LPBS(390,2),STDEVs_EPBS(390,3),STDEVs_LPBS(390,3); ...
              'CI95lower',STDEVs_EPBS(11,1),STDEVs_LPBS(11,1),STDEVs_EPBS(11,2),STDEVs_LPBS(11,2),STDEVs_EPBS(11,3),STDEVs_LPBS(11,3); ...
              'CI90upper',STDEVs_EPBS(380,1),STDEVs_LPBS(380,1),STDEVs_EPBS(380,2),STDEVs_LPBS(380,2),STDEVs_EPBS(380,3),STDEVs_LPBS(380,3); ...
              'CI90lower',STDEVs_EPBS(21,1),STDEVs_LPBS(21,1),STDEVs_EPBS(21,2),STDEVs_LPBS(21,2),STDEVs_EPBS(21,3),STDEVs_LPBS(21,3)}
RHOResults = {'','RHO(XBetaE,XBetaL)','RHO(etaE,etaL)','RHO(THETA_E,THETA_L)'; ...
              'PE',RHO(1),RHO(2),RHO(3); ...
              'StdErr',std(RHOPBS(:,1)),std(RHOPBS(:,2)),std(RHOPBS(:,3)); ...
              'CI95upper',RHOPBS(390,1),RHOPBS(390,2),RHOPBS(390,3); ...
              'CI95lower',RHOPBS(11,1),RHOPBS(11,2),RHOPBS(11,3); ...
              'CI90upper',RHOPBS(380,1),RHOPBS(380,2),RHOPBS(380,3); ...
              'CI90lower',RHOPBS(21,1),RHOPBS(21,2),RHOPBS(21,3); ...
              'Pval',2*normcdf(RHO(1)/std(RHOPBS(:,1))),2*normcdf(RHO(2)/std(RHOPBS(:,2))),2*normcdf(RHO(3)/std(RHOPBS(:,3)))} 
toc
clear betaregPBS STDEVs_EPBS1 STDEVs_EPBS2 STDEVs_EPBS3 STDEVs_LPBS1 STDEVs_LPBS2 STDEVs_LPBS3
clear RHOPBS1 RHOPBS2 RHOPBS3 STDEVs_EPBS STDEVs_LPBS RHOPBS STDEVs_E STDEVs_L RHO
clear NsimZ SimZ Cov_etaPBS Var_eta_EPBS Var_eta_LPBS XBeta_EPBS XBeta_LPBS Cov_XBetaPBS Var_XBetaPBS temp
clear BETA_EPBS BETA_LPBS muEPBS muLPBS sigmaEPBS sigmaE_fPBS sigmaE_blkPBS sigmaE_hspPBS SIGMAePBS 
clear sigmaLPBS sigmaL_fPBS sigmaL_blkPBS sigmaL_hspPBS SIGMAlPBS rhoELPBS rhoEL_fPBS rhoEL_blkPBS rhoEL_hspPBS RHOelPBS 



Y_Eimputese = nan(size(Y_Eimpute)) ;
Y_Limputese = nan(size(Y_Limpute)) ;
lY_Eimputese = nan(size(Y_Eimpute)) ;
lY_Limputese = nan(size(Y_Limpute)) ;
for ii=1:length(Y_Eimpute)
    Y_Eimputese(ii) = std(Y_EimputeBS(ii,:),'omitnan') ; 
    Y_Limputese(ii) = std(Y_LimputeBS(ii,:),'omitnan') ;
    lY_Eimputese(ii) = std(log(Y_EimputeBS(ii,:)),'omitnan') ;
    lY_Limputese(ii) = std(log(Y_LimputeBS(ii,:)),'omitnan') ;
end



format shortg

stdlThetaEShrunk  = sqrt( stdlThetaE^2 - mean(lY_Eimputese.^2) ) ;
ShrinkageEcompare = [stdlThetaE stdlThetaEShrunk 100*(stdlThetaEShrunk - stdlThetaE)/stdlThetaE] %#ok
ShrinkFactorE     = (stdlThetaEShrunk^2)./( stdlThetaEShrunk^2 + lY_Eimputese.^2 ) ;
Y_EimputeShrunk   = exp( log(Y_Eimpute).*ShrinkFactorE + mean(log(Y_Eimpute))*(1-ShrinkFactorE) ) ;

stdlThetaLShrunk  = sqrt( stdlThetaL^2 - mean(lY_Limputese.^2) ) ;
ShrinkageLcompare = [stdlThetaL stdlThetaLShrunk 100*(stdlThetaLShrunk - stdlThetaL)/stdlThetaL] %#ok
ShrinkFactorL     = (stdlThetaLShrunk^2)./( stdlThetaLShrunk^2 + lY_Limputese.^2 ) ;
Y_LimputeShrunk   = exp( log(Y_Limpute).*ShrinkFactorL + mean(log(Y_Limpute))*(1-ShrinkFactorL) ) ;



format long



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%            FIGURE 7:                      %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(1,2,1)
hold on
scatter(log(Y_Eimpute(Q>=2)),log(Y_EimputeShrunk(Q>=2)),'+')
scatter(log(Y_Eimpute(Q<2)),log(Y_EimputeShrunk(Q<2)),'^')
plot([-1.5,2.5],[-1.5,2.5])
box on
grid on
xlabel('RAW PRODUCTIVITY ESTIMATE $\log(\hat{\theta}_{p})$','Interpreter','latex','FontSize',20)
ylabel('$\log(\hat{\theta}_{p})$ EB SHRUNK PRODUCTIVITY ESTIMATE ','Interpreter','latex','FontSize',20)
legend('Active Students','Marginal/Inactive Students','45-Degree Line','FontSize',20,'FontWeight','bold','Location','northwest')
subplot(1,2,2)
hold on
scatter(log(Y_Limpute(Q>=2)),log(Y_LimputeShrunk(Q>=2)),'+')
scatter(log(Y_Limpute(Q<2)),log(Y_LimputeShrunk(Q<2)),'^')
plot([-10,2],[-10,2])
box on
grid on
xlabel('RAW MOTIVATION ESTIMATE $\log(\hat{\theta}_{m})$','Interpreter','latex','FontSize',20)
ylabel('$\log(\hat{\theta}_{m})$ EB SHRUNK MOTIVATION ESTIMATE ','Interpreter','latex','FontSize',20)
legend('Active Students','Marginal/Inactive Students','45-Degree Line','FontSize',20,'FontWeight','bold','Location','northwest')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%            FIGURE 4:                      %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
hold on
scatter(log(Y_EimputeShrunk(36*kappapre<11)),log(Y_LimputeShrunk(36*kappapre<11)),50,'vr')
scatter(log(Y_EimputeShrunk(36*kappapre>=11 & 36*kappapre<18)),log(Y_LimputeShrunk(36*kappapre>=11 & 36*kappapre<18)),50,'ob')
scatter(log(Y_EimputeShrunk(36*kappapre>=18)),log(Y_LimputeShrunk(36*kappapre>=18)),50,'pk')
scatter(mean(log(Y_EimputeShrunk(36*kappapre<11))),mean(log(Y_LimputeShrunk(36*kappapre<11))),300,'vr','filled')
scatter(mean(log(Y_EimputeShrunk(36*kappapre>=11 & 36*kappapre<18))),mean(log(Y_LimputeShrunk(36*kappapre>=11 & 36*kappapre<18))),300,'ob','filled')
scatter(mean(log(Y_EimputeShrunk(36*kappapre>=18))),mean(log(Y_LimputeShrunk(36*kappapre>=18))),450,'pk','filled')
legend('Pre-Exam Lower Tercile Students','Pre-Exam Middle Tercile Students','Pre-Exam Top Tercile Students','Location','NorthEast','FontWeight','bold','FontSize',20)
xlabel('log(\theta_p)','FontWeight','bold','FontSize',20)
ylabel('log(\theta_m)','FontWeight','bold','FontSize',20)
xlim([-1,2.2])
ylim([-6.5,0])
box on
grid on





if Model==4.05
    diary('.\STAGE4\InitialHCProduction_RF4point05.txt')
    Model  %#ok
    disp(datetime)
    run AAA_FILE05_InitialHCProductionRF.m
    diary off
elseif Model==4.085
    diary('.\STAGE4\InitialHCProduction_RF4point085.txt') 
    Model  %#ok
    disp(datetime)
    run AAA_FILE05_InitialHCProductionRF4point085.m
    diary off
end


if Model==4.05
    diary('.\STAGE4\HCProduction_ShortRun_RegOutput4point05.txt')
    Model  %#ok
    disp(datetime)
    run AAA_FILE06_HCProductionShortRun4point05.m
    diary off
elseif Model==4.085
    diary('.\STAGE4\HCProduction_ShortRun_RegOutput4point085.txt')
    Model  %#ok
    disp(datetime)
    run AAA_FILE06_HCProductionShortRun4point085.m
    diary off
end



